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ABSTRACT 

Much of the interstellar medium in disk galaxies is in the form of neutral 
atomic hydrogen, H I. This gas can be in thermal equilibrium at relatively low 
temperatures, T 300 K (the cold neutral medium, or CNM) or at temperatures 
somewhat less than 10 4 K (the warm neutral medium, or WNM). These two 
phases can coexist over a narrow range of pressures, P m i n < P < P max - We 
determine P m i n and P max in the plane of the Galaxy as a function of Galactocentric 
radius R using recent determinations of the gas heating rate and the gas phase 
abundances of interstellar gas. We provide an analytic approximation for P m j n 
as a function of metallicity, far-ultraviolet radiation field, and the ionization rate 
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of atomic hydrogen. Our analytic results show that the existence of P m i n , or 
the possibility of a two-phase equilibrium, generally requires that H + exceed 
C + in abundance at P m i n . The abundance of H + is set by EUV/soft X-ray 
photoionization and by recombination with negatively charged PAHs. In order 
to assess whether thermal or pressure equilibrium is a realistic assumption, we 
define a parameter T = t coo iAshk where t coo \ is the gas cooling time and t s hk is 
the characteristic shock time or "time between shocks in a turbulent medium" . 
For T < 1 gas has time to reach thermal balance between supernovae induced 
shocks. We find that this condition is satisfied in the Galactic disk, and thus 
the two-phase description of the interstellar H I is approximately valid even in 
the presence of interstellar turbulence. Observationally, the mean density (nni) 
is often better determined than the local density, and we cast our results in 
terms of (nni) as well. Over most of the disk of the Galaxy, the H I must be in 
two phases: the weight of the H I in the gravitational potential of the Galaxy 
is large enough to generate thermal pressures exceeding P m i n , so that turbulent 
pressure fluctuations can produce cold gas that is thermally stable; and the mean 
density of the H I is too low for the gas to be all CNM. Our models predict the 
presence of CNM gas to R ~ 16 — 18 kpc, somewhat farther than previous 
estimates. An estimate for the typical thermal pressure in the Galactic plane for 
3 kpc ^ R Z 18 kpc is P th /k ~ 1.4 x 10 4 exp(-P/5.5 kpc) K cm 3 . At the solar 
circle, this gives Pth/k ~ 3000 K cm -3 . We show that this pressure is consistent 
with the C I*/C I to t ratio observed by Jenkins & Tripp (2001) and the CNM 
temperature found by Heiles & Troland (2002). 

We also examine the potential impact of turbulent heating on our results and 
provide parameterized expressions for the heating rate as a function of Galactic 
radius. Although the uncertainties are large, our models predict that including 
turbulent heating does not significantly change our results and that thermal 
pressures remain above P m i n to R ~ 18 kpc. 

Subject headings: ISM: clouds — ISM: general — ISM: structure 



1. INTRODUCTION 

The interstellar medium (ISM) has a complex thermal and ionization structure. Much of 
the neutral atomic gas is observed to be either warm neutral medium (WNM) with T ~ 10 4 K 
or cold neutral medium (CNM) with T ~ 100 K (Kulkarni & Heiles 1987; Dickey & Lockman 
1990). Some of the warm gas is partially ionized, the warm ionized medium (WIM), which 
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also has T ~ 10 4 K (McKee & Ostriker 1977; Reynolds 1983; Haffner, Reynolds, & Tufte 
1999). A small mass fraction of the gas is in the form of hot ionized medium (HIM) with 
T ~ 10 6 K (Cox k Smith 1974; McKee & Ostriker 1977). Inside the solar circle, about half 
the interstellar gas is molecular (Scoville & Sanders 1987; Bronfman et al. 1988, 2000). 

A significant simplification occurs if one focuses on the neutral atomic gas, the CNM 
and WNM. Some decades ago, Field, Goldsmith, & Habing (1969) demonstrated that the 
CNM and WNM could coexist in pressure equilibrium, so that the neutral atomic gas could 
be considered to be a two-phase medium. They assumed that cosmic rays dominate the 
heating, but it was subsequently realized that UV starlight dominates the heating due to 
photoelectric emission from the dust grains in the gas (Watson 1972). Using the photoelectric 
heating rates calculated by Bakes & Tielens (1994), Wolfire et al. (1995; hereafter WHMTB) 
investigated the thermal balance of the WNM and CNM phases in the local ISM and showed 
that the two-phase model is in good agreement with a wide variety of data on the ISM in 
the solar vicinity. 

What is the evidence for a two-phase medium elsewhere in the Galaxy? In the inner 
Galaxy, Garwood & Dickey (1989) found that there is H I emission (which can originate from 
both CNM and WNM) at all velocities allowed by Galactic rotation. On the other hand, they 
found that absorption (which originates only from CNM at their sensitivity) is somewhat 
less pervasive, particularly within 2 kpc of the Galactic Center. Liszt, Braun, & Greisen 
(1993) suggested that the H I absorption in the inner Galaxy at R > 2 kpc is much higher 
than that reported by Garwood & Dickey (1989). Kolpak et al. (2002) recently repeated 
the earlier H I absorption study and confirmed the presence of cold gas in the inner Galaxy 
with an absorption coefficient at R = 5 kpc approximately 5 times higher than reported 
by Garwood & Dickey (1989). In the outer Galaxy, the presence of WNM is reasonably 
well established (Kulkarni & Heiles 1987), whereas that of a widely distributed CNM is less 
so. Carilli, Dwarakanath, & Goss (1998) have measured the temperature of the WNM in 
absorption features seen towards Cygnus A at distances of 9 kpc and 12 kpc (and z height of 
~ 1 kpc) using the Westerbork radio telescope. They find gas temperatures of ~ 6000 K and 
~ 4800 K respectively, which are consistent with the low pressure and low UV field models 
of WHMTB for atomic gas above the plane. Several high-velocity absorption components 
have been observed in H I (Colgan, Salpeter, & Terzian 1988) and Na I (Sembach & Danks 
1994) that arise from CNM clouds at Galactic radii R < 14 kpc. Kolpak et al. (2002) show 
H I absorption to R ~ 17 kpc. Molecular clouds, which presumably form from the CNM 
phase, are traced to at least R ~ 20 kpc (Wouterloot & Brand 1989; Heyer et al. 2001) 
with an extremely distant H II region and molecular cloud complex at R = 28 kpc (Digel, de 
Geus, & Thaddeus 1994). Wouterloot et al. (1990) and Bronfman et al. (2000) show that the 
molecular surface density can be fit by a radial exponential in the outer Galaxy to ~ 18 kpc 
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(see also Williams & McKee 1997). In their study of the Perseus arm, Heyer et al. (1998), 
however, find that the molecular gas disk is effectively truncated at R ~ 13.5 kpc. These 
two results could be consistent if molecular gas extends to greater radii in directions other 
than those studied by Heyer et al. (1998) or if isolated molecular clouds extend to distances 
much greater than the molecular surface density can be reliably measured from CO surveys. 

Thus direct observations of cold H I or molecular gas extend to at least R ~ 18 kpc. Star 
formation provides an indirect test for the presence of CNM, since the gas that forms stars 
presumably goes through the stage of being cold H I. In the Galaxy, near-infrared sources, 
IRAS sources, and H II regions (Wouterloot et al. 1990; Rudolph et al. 1996; Kobayashi & 
Tokunaga 2000; Snell, Carpenter, & Heyer 2002) are seen out to R ~ 17 — 20 kpc, suggesting 
that CNM extends out to at least that distance. 

What can be learned from observations of H I in other galaxies? Two-phase atomic 
gas has been observed using 21 cm absorption techniques in several extragalactic systems 
including M31 (Dickey & Brinks 1993; Braun & Walterbos 1992), M33 (Dickey & Brinks 
1993), the LMC (Mebold et al. 1997), and the SMC (Dickey et al. 2000). Braun (1997), using 
the VLA, examined the neutral hydrogen emission in 11 nearby spirals. By associating high 
brightness, narrow emission components with cold gas, he finds that the fraction of cold gas 
remains relatively constant until the B band surface brightness falls to fiB ~ 25 mag arcsec -2 , 
i.e., the R25 radius. At larger radii, the fraction drops, although in some systems more than 
10% of the H I is in the form of cold gas out to (1.5 — 2) x R 2 ^. Since the extinction-corrected 
radius of the Galaxy is R25 = 12.25 kpc (de Vaucouleurs 1983), cold gas in the Galaxy could 
extend out to 1.5 - 2 x R 2 n, or R ~ 18.4 - 24.5 kpc. Sellwood & Balbus (1999) interpret the 
line width of H I in the outer parts of galaxies as being due to CNM that is being stirred 
by the magnetorotational instability. Evidence for recent star formation in the outer disk of 
M31 is presented by Cuillandre et al. (2001), who find a population of B stars at 1.7 x i? 2 5, 
pointing towards the presence of cold gas in the outer parts of galaxies. 

The fact that the neutral atomic gas in the Galactic ISM is in two phases is a powerful 
result, since two phases can coexist only over a narrow range of pressure, P max > P > -P m in 
with P max ^ 3P min (Field, Goldsmith, & Habing 1969, see § 6). It is thus possible to estimate 
the thermal pressure of the H I with reasonable accuracy — when it is in two phases — from 
knowing the gas phase abundances, the dust properties, and the intensity of the radiation 
field. We used this property to study gas in the Galactic halo and constrain the properties of 
the High Velocity Clouds in Wolfire et al. (1995b). Corbelli & Salpeter (1988) have argued 
that achieving the condition for a two-phase equilibrium is a necessary step in initiating 
star formation in young galaxies, while Parravano (1988, 1989) has suggested that two-phase 
equilibria play a key role in regulating the rate of star formation in disk galaxies. 
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The primary goal of this paper is to predict the average thermal pressure of the ISM 
as a function of position in the Galaxy using the two-phase criteria. To do this, we shall 
extend the models of WHMTB to the inner and outer Galaxy. In light of the observational 
evidence that cold gas exists in the outer Galaxy, we shall carry out our model calculations 
for Galactic radii between 3 kpc and 18 kpc. Knowing the thermal pressure allows one to 
predict the intensities of the dominant cooling lines of the gas, such as that of C II 158 /im, 
and examine the heating and cooling routes which determine the energy budget. Locally, the 
thermal pressure in the ISM is measured through ultraviolet absorption line studies (Jenkins 
et al. 1983; Jenkins & Tripp 2001). In the near future, telescopes such as ALMA, SOFIA, 
and Herschel will provide additional measurements of the thermal pressure and dominant 
cooling lines throughout the Galaxy and in other galaxies. These will test our model for the 
gas thermal balance and check the importance of thermal instability. 

We also can calculate whether the ISM could exist as pure WNM at various positions 
in the Galaxy by comparing the weight of the H I layer with P max . The problem of deter- 
mining the phase structure of the H I in the outer Galaxy has been considered previously by 
Elmegreen & Parravano (1994), who find a transition to mainly WNM at R ~ 12 — 14 kpc. 
Our results are compared with theirs in § 8.3. 

Although the focus of this paper is on the determination of the thermal pressure in 
Galactic H I, it is well known that the thermal pressure is only a small part of the total 
pressure in the gas; in particular, the turbulent pressure is considerably greater than the 
thermal pressure (Boulares & Cox 1990). In § 2, we discuss the relation between the turbulent 
pressure and the thermal pressure and determine the conditions under which it makes sense 
to consider multi-phase equilibria in a turbulent medium. We also discuss in Appendix B 
the dissipation of turbulent energy in the ISM and its potential effects on our results. In 
§ 3 we discuss the distribution of gas and dust in the Galaxy between these radii, together 
with the abundances we have adopted. The heating and ionization in the gas are governed 
by energetic photons and particles, which are discussed in § 4. The thermal and chemical 
processes in our model are slightly modified from those discussed by WHMTB; the differences 
are briefly described in § 5. The results of our calculations are presented in § 6. We then 
construct a simple analytic model of a two-phase equilibrium that shows how the properties 
of the equilibrium scale with the input parameters (§ 7). We compare our model with local 
and extragalactic observations in § 8 and discuss our results in § 9. 
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2. Turbulent Pressure in the Multiphase Interstellar Medium 

The ISM is observed to be highly turbulent (Larson 1979). When averaged over the 
vertical structure of the galactic disk, the turbulent velocity dispersion exceeds the thermal 
velocity dispersion, and correspondingly the turbulent pressure exceeds the thermal pressure 
(Boulares & Cox 1990). Theoretical arguments (Spitzer 1968, 1978; McKee & Ostriker 1977) 
and numerical simulations (Bregman, Parriott, & Rosen 1999; Korpi et al. 1999) show that 
these turbulent motions can be accounted for by the injection of energy by supernovae. 

Insofar as the disk of the Galaxy is in approximate hydrostatic equilibrium, the total 
pressure in the midplane of the disk must balance the weight of the material above it (Parker 
1969; Boulares & Cox 1990), and the turbulent motions are an important contributor to 
the total pressure. In the solar neighborhood, the total pressure at the midplane is about 
P/k ~ 2.8 x 10 4 K cm -3 (Boulares & Cox 1990), about 10 times greater than the median 
thermal pressure of ~ 3000 K cm -3 (Jenkins & Tripp 2001). 1 A similar situation occurs in 
molecular clouds, where the total pressure can also be about an order of magnitude greater 
than the thermal pressure (depending on scale — see Larson 1981). If the thermal pressure is 
such a small fraction of the total, why are arguments based on thermal pressure equilibrium 
relevant in the interstellar medium? 

The answer is that the turbulent motions determine the temporal and spatial structure of 
the thermal pressure. A key feature of turbulence is that the motions are spatially correlated, 
so that the rms velocity increases with scale — i.e., there is a relation between the line width 
and the size of a region, 

O = <7(1K PC , (1) 

where cr(l) is the ID turbulent velocity on the scale of 1 pc and £ pc is the scale in units of 
parsecs. Larson's (1979) data on H I clouds give <r(l) = 0.64 km s _1 and q = 0.37 over the 
range 1 ~ £ pc ~ 1000. (Larson 1981 subsequently found a similar line width-size relation 
for molecular clouds.) We can estimate cr(l) from a more homogeneous data set by using 
the recent H I study by Heiles & Troland (2002). They find that the mass- weighted velocity 
dispersion is 7.1 km s _1 for CNM clouds and 11.4 km s _1 for WNM clouds. The typical line 
of sight in this survey is at a Galactic latitude slightly greater than 30°, so the typical path 
length is about twice the half- height of the disk — i.e., the full disk thickness. From Dickey & 
Lockman's (1990) model for the vertical distribution of the H I, we infer that the FWHM of 
the CNM is 212 pc, whereas that for the WNM is 530 pc. We assume that the outer scale of 



1 The mean pressure quoted by Jenkins & Tripp (2001) is P/k — 2240 K cm -3 based on data from the 
Space Telescope Imaging Spectrograph. In § 8 we find that corrections for gas temperature and atomic 
constants raise the pressure derived from the observations to <~ 3000 K cm~ 3 . 
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the turbulence is greater than the disk thickness, so that the power-law behavior extends over 
this length scale; this assumption is consistent with the results of Larson (1979), who did not 
find a break in the power-law behavior out to 1 kpc in the Galaxy, and with the results of 
Lazarian & Pogosyan (2000), who did not find a break between 40 pc and 4 kpc in the data 
for the Small Magellanic Cloud. We shall set q = 1/3, as expected for subsonic turbulence. 
On scales large enough for the turbulence to be supersonic, q is expected to be somewhat 
larger: Larson (1979) suggested that q would approach |, whereas Boldyrev (2002) finds 
q ~ 0.37 — 0.38 (where a is interpreted as an rms velocity). For the CNM, we then find 
er(l) = 1.2 km s -1 , whereas for the WNM we find <r(l) = 1.4 km s' 1 . The value of q is likely 
to be closer to the assumed 1/3 on the smaller length scales associated with the CNM, so 
we shall adopt a(l) = 1.2 km s _1 for the H I in the solar neighborhood. Although this value 
was determined from the velocity dispersion of different H I clouds, we shall assume that 
it applies within individual H I clouds as well. We note that this assumption is consistent 
with the results of Linsky & Wood (1996), who found that the velocity dispersion within the 
H I toward the nearby star a Cen is about 1.2 km s" 1 , corresponding to er(l) = 1.1 km s _1 . 
The use of a single turbulent velocity law in both the CNM and the WNM is undoubtedly 
an oversimplification; in particular, there is a range of scales over which the turbulence is 
supersonic in the CNM and subsonic in the WNM, and q might be expected to have different 
values within individual clouds of CNM and WNM over this range. 

Let £p be the scale on which the turbulent pressure begins to dominate the thermal 
pressure; on scales less than £ P the gas typically will be in approximate thermal pressure 
equilibrium. To estimate £p, we equate the thermal velocity dispersion, cr th = 0.80T 2 1//2 km 
s _1 , with the turbulent velocity dispersion, so that 



o~th 

(7(1) 



pc - 0.3T 2 3/2 pc, (2) 



where T 2 = T/(100 K) and where the numerical evaluation is for our fiducial case. Because 
lp depends on the cube of the uncertain quantity er(l), the numerical value of £p is quite 
uncertain. Bearing this in mind, we find that the CNM should be in approximate thermal 
pressure equilibrium on scales £ < £ P (CNM)~ 0.3 pc. For the WNM, with a typical tem- 
perature of about 8000 K, we find f P (WNM) ~ 215 pc. This is somewhat larger than the 
size of the turbulent cells in the warm gas, 60 pc, found in the 3D numerical simulations of 
the ISM by Korpi et al. (1999). 

In a multiphase medium, the CNM is embedded in warm gas, either predominantly 
neutral (WNM) or ionized (WIM), which in turn may be embedded in hot gas, the HIM. 
Since the sound speed in the WIM exceeds that in the WNM, £ P (WIM) is somewhat larger 
than £p(WNM); both are much larger than the typical size of a CNM cloud, which is ~ 1 — 2 
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pc (McKee & Ostriker 1977). As a result, the CNM clouds should typically be embedded in 
a medium that is spatially isobaric. The surface layers of the cloud [i.e., those layers within 
a distance £p(CNM) of the surface] should typically have the same thermal pressure as the 
ambient warm medium, whereas the inner regions of the cloud should have the time-averaged 
value of the ambient thermal pressure. Elmegreen (1997) shows that in a ID simulation of 
interacting magnetized clouds the gas maintains approximate phase equilibrium with the 
intercloud medium. 

This description of the relative roles of thermal and turbulent pressures in the ISM is 
consistent with the three-phase model of the ISM by McKee & Ostriker (1977). In their 
model, the supernova-generated HIM is pervasive and sets the thermal pressure in the em- 
bedded clouds of CNM and WNM. The ISM is viewed as the superposition of many supernova 
remnants; as a result there are substantial fluctuations in the thermal pressure, in qualita- 
tive agreement with observations (Jenkins et al. 1983; Jenkins & Tripp 2001). The turbulent 
motions provided by supernovae are consistent with observations of H I velocities in the 
Galactic plane. Turbulent pressure is the single largest contributor to the support of the 
ISM in the gravitational field of the disk (Boulares & Cox 1990), and McKee (1990) has 
argued that the turbulent motions that produce this pressure are produced by supernovae. 
Slavin & Cox (1993) and Cox (1995) have argued that magnetic fields in the ISM limit the 
size of supernova remnants so that they occupy only a small fraction of the volume. Such a 
model cannot account for the observed level of turbulence, however. 

Under what conditions is the turbulence in the ISM weak enough that a two-phase 
description of the H I is valid? We first consider the cooling time t coo i for an element of gas 
to return to thermal equilibrium after a significant perturbation. From our numerical results 
presented in § 6, we find the cooling rate at a constant thermal pressure of Pth/k = 3000 K 
cm" 3 is given by 

A = 5.7 x 10" 26 T 4 a8 erg cm 3 s~\ (3) 

where T 4 = T/(10 4 K) and where the fit is accurate to within a factor 1.35 for temperatures 
between T = 55 K and T = 8500 K. With this cooling rate, t coo i is given by 



f-cool 



\l.lnkT\ Rio/ Pth/k \ 

— ).-.-xl0^( 3000 ^ m _ 3 ) y, (4) 



-0.8 



where n is the hydrogen nucleus density (n = n H i + n H + + 2nn 2 )- The fit to the pressure 
term results in an overall accuracy within a factor 1.5 for temperatures between T = 55 K 
and T = 8000 K and thermal pressures between P th /k = 1000 K cm" 3 and P th /k = 6000 
K cm" 3 . In order for a two-phase description to be valid, this time must be less than 
the characteristic time t s hk for a shock to induce motions of order a t h or, equivalently (for 
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an isothermal shock), to double the pressure. 2 One method of estimating t s hk is based on 
observations of interstellar turbulence. In order to induce motions of a t h, the velocity of an 
isothermal shock must be y/2a t \ 1 . The time interval between shocks in a turbulent medium is 
then t s hk ~ ^p/\/20"th) where, as discussed above, the value of tp is quite uncertain. Defining 
a dimensionless parameter T that measures the strength of the turbulence, we obtain 



Y _ ^cool ^ q 2^0.2 / Pth/k 



tshk ' 4 VSOOOKcm- 3 



-0.8 r 

m ; i 

(5) 



(7(1) 



1.2 km s" 1 



So long as T ^ 1, a two-phase description of the H I is approximately valid, since the gas 
has time to reach thermal balance between shocks. This condition is well satisfied for the 
CNM (T ~ 0.1), but only marginally satisfied for the WNM (T ~ 0.3). 

We can also estimate t s hk and T analytically under the assumption that the shocks 
are produced by supernova remnants, although this estimate is necessarily uncertain. Let 
M shk (> v s ) be the mass of interstellar gas per supernova that suffers a shock with a shock 
velocity greater than v s . Since the shock must have a velocity \^2a t h in order to induce 
motions of <7 t h, we have 

^ ^ ^WNM ^ 

?SN,efT-^shk(> V^20"th) 

where Ewnm is the surface density of WNM gas, and ^sn, eff is the effective supernova rate per 
unit area. Allowing for the clustering of supernovae in associations, McKee (1989) estimated 
that the effective supernova rate is only about 36% of the actual rate. For a two-phase ISM, 
he also estimated M shk = 2460n^° N 1 M (t> s /100 km s -1 ) -9 / 7 M Q , where «wnm is the H nucleus 
density in the WNM and where we have assumed that a typical supernova remnant has an 
energy of 10 51 erg. (The results for a three-phase ISM with a substantial WNM filling factor 
should not differ qualitatively from the two-phase results.) At the solar circle, the supernova 
rate per unit area is about 3.8 x 10~ n SN pc~ 2 yr _1 (McKee & Williams 1997), so this 
gives t s hk — 5.3 x 10 6 yr for a density t^wnm — 0.3 cm~ 3 . The corresponding value of the 
turbulence parameter for the WNM is T ~ 0.9, somewhat larger than the value estimated 
above (and indicative of the uncertainty in this parameter). Both estimates of T suggest 
that the WNM will often exhibit non-equilibrium temperatures. Heiles & Troland (2002) 
find that about half the WNM is at temperatures below the minimum equilibrium value, 
which is qualitatively consistent with our estimate of T. 



2 In Siting equations (3) and (4) we have assumed ionization equilibrium holds although the recombination 
time in WNM gas is found to be comparable to, or up to 2 times greater than the cooling time. Since the 
cooling time in the WNM is proportional to l/n e for Lya cooling, and proportional to 6 for electrons 
recombining on to positively charged grains (see § 5), any lingering ionization after the passage of a shock 
will tend to decrease the cooling time. Thus, our fits provide a maximum cooling time in WNM gas. 
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Recent papers by Mac Low et al. (2001), Vazquez-Semadeni, Gazol, & Scalo (2000), 
and Vazquez-Semadeni et al. (2002) argue that thermal instability is of less importance 
in determining the distribution of phases than is turbulent dynamics followed by cooling 
towards thermal equilibrium. To some extent, the difference between these conclusions 
and those in the present paper are a matter of perspective: these authors focus on the 
time-dependent aspects of a turbulent medium, whereas we emphasize the utility of the 
equilibrium aspects. (We also note that Mac Low et al. (2001) did not include heating sources 
in their calculations.) In terms of the turbulence parameter T introduced above, we choose 
to approximate the case T ~ 1 for the WNM with the two-phase T 1 results, whereas they 
prefer to emphasize the case T»l. Also, in contrast to Mac Low et al., who concentrate 
on evaluating the thermal pressure distribution in the solar neighborhood, the focus of our 
work is on determining the mean (volume-averaged) thermal pressure in the ISM throughout 
the Galaxy, P t h,avc- Unless P t h, avc is in, or close to, the range P min ^ Pth.avc ~ Pmax, the gas 
will be almost all CNM or all WNM, and the complex interplay between cold and warm gas 
seen by these authors will not occur. By determining P min and P max , we can determine the 
range of pressures in which cold and warm gas can coexist; furthermore, we shall argue that 
most of the Galactic disk must have a thermal pressure such that this is the case. 



3. Gas and Dust in the Milky Way 

3.1. Distribution of H I 

We require the azimuthally averaged H I surface density distribution, E H i(P), the half 
width to half maximum of the H I emission above the plane, H^ J (R), and the mean H I 
density in the Galactic plane, (nni(P)) oc Shi(P)/P^ : (P). These will be used to trace 
the opacity in order to determine the distribution of energetic photons and particles in the 
Galaxy as discussed in § 4. 

Several surveys have been conducted of the Galactic H I distribution (e.g., Weaver & 
Williams 1973; Burton 1985; Kerr et al. 1986; Stark et al. 1992; Hartmann & Burton 1997). 
Lockman (2002) discusses several limitations in interpreting the observations including: ve- 
locity crowding which renders Shi(P) highly sensitive to the adopted rotation curve, optical 
depth effects which introduce uncertainty in the volume and surface densities, and a low 
dynamic range in emission, which makes small variations in the brightness temperature dif- 
ficult to measure. Thus, the distribution in Ehi(P) is not well determined. In addition, the 
derived scale height H^ 1 may depend on the method of analysis which could be preferen- 
tially sensitive to the CNM or WNM component. With these caveats in mind, we use the 
published H I data to derive our distributions. 
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The distribution of H I in the outer Galaxy was presented by Wouterloot et al. (1990) 
who combined the northern hemisphere (Kerr et al. 1986; Burton 1985) and southern hemi- 
sphere (Kerr et al. 1986) data. They provide plots of the H I surface density, Ehi(-R), from 
the 2nd and 3rd Galactic quadrants (northern and southern data) as well as the average 
surface density and find that the average radial distribution at R > 13 kpc is well fit by an 
exponential, Shi(-R) oc exp(— R/H^ 1 ), with H^ 1 = 4 kpc. A relatively flat rotation curve 
was used, similar to that derived by Brand & Blitz (1993). 

Dame (1993) and Lockman (1988) showed that minor (~ 2%) differences in the rotation 
curve can greatly affect the surface density derived from the data with differences in the 
surface density amounting to ~ 50%. Based on the northern hemisphere data, Dame (1993) 
found a peak in surface density near 12 kpc for a flat rotation curve, but a nearly constant 
surface density out to 17 kpc when using the slightly rising rotation curve of Kulkarni, Blitz, 
& Heiles (1982). Wouterloot et al. (1990) attributed the peak in the H I surface density to 
gas associated with the Perseus arm and show that the average of the northern and southern 
data partially smoothes the 12 kpc peak. 

In addition to the surface density enhancement at 12 kpc, the Wouterloot et al. (1990) 
plots of the outer Galaxy show a dramatic rise as R decreases towards the solar circle, 
reaching Ehi = 8.6 M & pc -2 at R = 9 kpc. (Note that in our notation Ehi does not include 
the mass associated with helium or the metals.) This surface density is inconsistent with the 
value of Shi(-Ro) — 5 M pc~ 2 at the solar circle determined by Dickey & Lockman (1990). 
Furthermore, the average surface density in the outer Galaxy does not seem to join smoothly 
onto the inner Galaxy where Shi is also 5 M Q pc~ 2 . The inferred pile up of H I towards 
the edge of the solar circle may be an artifact of a flat rotation curve with strictly circular 
velocities, and non-circular motions may alleviate this problem (Blitz & Spergel 1991). An 
example of the H I surface density retaining the peak at R is shown in Oiling & Merrifield 
(1998), who also argue for R ~ 7.1 kpc. In our post-Copernican world, we believe that 
a narrow density enhancement centered at R = 8.5 kpc is unrealistic, and have smoothly 
joined the Wouterloot et al. (1990) surface densities at R > 10 kpc onto the distribution 
for R < 8.5 kpc. We are mainly concerned with the disk properties in the inner and outer 
Galaxy, and thus errors less than a factor 2 in a thin region between 8.5 and 10 kpc are not 
critical. 

The observational evidence suggests that the H I surface density is constant in the inner 
Galaxy between 8.5 kpc and ~ 4 kpc and then drops by a factor ~ 3 by 1.5 kpc (Dickey & 
Lockman 1990; Liszt 1992; Dame 1993). Our piecewise analytic fit to the H I surface density 
data is shown in Figure 1 and is given by 



EhiCR) = 1.4/4-0.6 M pc~ 2 



(3 < R k < 4), 



(7) 
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= 5 M pcT 2 

= -1.12 + [6.12 x (74/8-5)] M pc 
= 8.24 e -(^- 13 )/ 4 M pc- 2 



-2 



(4 < < 8.5), 
(8.5 <R k < 13), 
(13 < i? fe < 24), 



(8) 
(9) 
(10) 



where Rk = i2/(l kpc) and where we use the conversion 1 M pc -2 = 1.25 x 10 20 H I cm -2 . 
Dame (1993) found that using the H I surface density given in Liszt (1992), the H I mass 
at R < R (excluding the Galactic center) is M m ~ 1.7 x 10 9 M Q while our distribution 
gives a mass M H i ~ 1.0 x 10 9 M in this same range. A simple estimate of the maximum 
H I mass is M H i < 7ri2 2 E H i(i2o) ~ 1.2 x 10 9 M and our lower estimate seems reasonable. 
We find an H I mass of Mhi ~ 4.3 x 10 9 M in the range 8.5 < R < 18 kpc consistent with 
Wouterloot et al. (1990) who report M HI = 5.3 x 10 9 M between R = 8.5 kpc and 24 kpc. 

The H I half height, Hf 1 (R), at R > R is taken from the data of Wouterloot et al. 
(1990). As for the surface density distribution, we assume that the data at 10 kpc should 
be smoothly joined to that at 8.5 kpc. With this assumption we find that the variation in 
the height over the entire range from 8.5 kpc to 18 kpc can be reasonably well fit by an 
exponential of scale length ~ 6.7 kpc. From Dickey & Lockman (1990) we assume the half 
height at R < R is approximately constant and given by Hf l (R < R ) = 115 pc. We 
note that Malhotra (1995) finds a constant half height of ~ 118 pc in the inner (R < 5.1 
kpc) Galaxy which rises to ~ 260 pc at the solar circle; a height which is nearly identical 
to that of the WNM (~ 265 pc) found by Dickey & Lockman (1990). Note that Malhotra 
(1995) fit both the midplane height and the height of the H I above the midplane and thus 
possible affects due to disk corrugation are removed from her results. The dramatic rise 
in height at the solar circle might be partly attributed to the correction for saturated H I 
emission by Dickey & Lockman (1990) which was not accounted for by Malhotra. Including 
this correction increases the intensity of the CNM component and tends to weight the height 
more towards the CNM height than the WNM height. 

Figure 2 shows the H I height distribution with an analytic fit given by 



The H I surface density and height are used to scale the mean midplane density (nni(-R)) oc 
^hi(-R) / H^ l {R). Liszt (1992) notes that with a density of (nni) ~ 0.4 cm -3 and a single 
temperature T spin = 135 K, both the emission and absorption characteristics of the H I in the 
plane can be simultaneously modeled. On the other hand, he points out that a higher den- 
sity is required to match the observed H I surface density of Dickey & Lockman (1990) who 
find (nni) — 0.57 cm -3 . With the mean density of hydrogen nuclei as derived from extinc- 
tion studies n H i + 2n H2 = 1.15 H-nuclei cm -3 (Bohlin, Savage, & Drake 1978), the Dickey & 



Hf\R) = 115 pc 

= 115 e ( R ^8.5)/6.7 



(3 < R k < 8.5), 
pc (8.5 < R k < 18). 



(11) 
(12) 
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Lockman (1990) mean H I density requires a molecular density of 2nn 2 ~ 0.58 cm -3 , roughly 
consistent with the observations of Bronfman et al. (1988) who find 2nn 2 ~ 0.5 cm -3 . We 
adopt a local mean H I density of (n H i) = 0.57 cm -3 and show the run of mean midplane 
H I density with R in Figure 3. 

Our H I distribution differs from that adopted by Ferriere (1998) for her global models 
of the Galactic ISM. Ferriere adopted a constant surface density Shi(-R) = 5 M Q pc~ 2 to 
R = 20 kpc with a scale height that increases linearly with R at R > 8.5 kpc and a mean 
midplane density that decreases as 1/R. Our surface density is thus equal to, or higher than 
Ferriere (1998) out to a radius of R = 15 kpc and then drops to lower values. Our midplane 
density is slightly higher until R = 14 kpc and then drops to a factor of ~ 4 lower at R = 18 
kpc. 



3.2. Distribution of H 2 

The distribution of the H2 surface density in the Galaxy is required primarily as a cosmic 
ray opacity source. We shall also use the distribution of H2 height to estimate the height of 
OB stars in the Galactic plane. We use a surface density distribution from Bronfman et al. 
(2000) that is a Gaussian in Galactic radius in the inner Galaxy and a radial exponential in 
the outer Galaxy. The distribution has a peak value of 4.5 M pc~ 2 centered on R = 4.85 
kpc, a full width at half maximum extent equal to 4.42 kpc, and a radial exponential of 
scale length H^ 2 = 2.89 kpc beyond R > 6.97 kpc. The H 2 surface density at the solar 
circle is £h 2 (-Ro) — 1-4 M Q pc~ 2 . These values are derived from an average of the data 
from both northern and southern Galactic quadrants. Bronfman et al. (2000) find an H 2 
mass of M(H 2 ) = 6.1 x 10 8 M for 1.7 < R < 8.5 kpc, and M(H 2 ) = 2 x 10 8 M Q for 
8.5 < R < 14.5 kpc. By comparison, Williams & McKee (1997) adopted a simpler, purely 
exponential distribution for the H 2 and found M(H 2 ) = 7.1 x 10 8 M for R between 1.7 kpc 
and 8.5 kpc (their quoted value of 1.0 x 10 9 M Q included the He mass). Figure 1 shows the 
H 2 surface density distribution. 

In the inner Galaxy, Bronfman et al. (2000) find the H 2 height to half maximum density 
is roughly constant at Hf 2 (R < R ) ~ 59 pc. In the outer Galaxy we simply scale the H 2 
height by the H I height of Wouterloot et al. (1990) (eq. [12]). 
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3.3. Distribution of Dust and Metals 

The photoelectric heating rate depends on the abundance of the smallest dust particles, 
which may primarily consist of polycyclic aromatic hydrocarbons or "PAHs", while the 
radiative cooling rate depends mainly on the gas phase abundances of carbon and oxygen. 
The photoelectric heating rate has been studied in detail by Bakes & Tielens (1994) using 
a realistic model for the PAH and grain size distribution. For the grains (i.e., particles 
with a > 15 A), they adopted the MRN grain size distribution for spherical grains (i.e., 
n[a]da oc a~ 3 ' 5 da; Mathis et al. 1977). The PAH molecules (a < 15 A) are assumed to be 
small disks. Their distribution is given by n(N c )dN c = 1.15 x 10 -5 iVc 2 - 25 djVc (with Ac the 
number of C-atoms in the molecule), which ensures that the total volume in disks between 12 
and 275 C-atoms is equal to that of the MRN spherical grain size distribution between 3 and 
15 A. WHMTB have shown that this heating rate can adequately heat the local WNM and 
CNM phases. The characteristics of the PAH distribution adopted by WHMTB correspond 
to a total fraction of C in the form of PAHs of 14 x 10~ 6 relative to H. Recent analysis of 
ISO observations of the PAH emission in the galaxy conclude that this fraction is actually 
22 x 10" 6 relative to H (Tielens et al. 1999). We note that this corresponds to a total PAH 
abundance, Tip ah/ ^ of 6 x 1CT 7 as compared to 4 x 1CT 7 in WHMTB. We have adopted 
the former, higher value, which enters in the photoelectric heating rate and the chemical 
network. As will be discussed in § 8, these new values still reproduce the observed [C II] 158 
fim cooling rates and the thermal pressures in the local ISM. 

We use a local gas phase abundance of metals consistent with the Hubble Space Telescope 
observations of Sofia et al. (1997), Cardelli et al. (1996), and Meyer, Jura, & Cardelli (1998) 
with ric/n = 1.4 x 10~ 4 for carbon and rio/n = 3.2 x 10~ 4 for oxygen. The gas phase 
abundances of these elements are seen to be independent of the physical conditions in the 
diffuse medium, showing relatively constant values for fractional H 2 abundances between 
lognH 2 /[^H 2 + ^hi] = —5.0 and —0.2. Although of relatively minor importance in these 
calculations, we also include Si, S, Mg, and Fe and, with the exception of S, increase their 
depletions with gas density as in Jenkins (1987) and Van Steenberg & Shull (1988) (see also 
WHMTB). 

The optical and infrared line diagnostics seem to be converging on a value for the 
oxygen abundance gradient in the Galaxy that is roughly consistent with the early findings 
of Shaver et al. (1983), which were based on radio observations of the temperatures of H II 
regions (e.g., Torres- Peimbert & Peimbert 1977; Simpson et al. 1995; Afnerbach et al. 1996, 
1997; Gummersbach et al. 1998; Rolleston et al. 2000). We will use oxygen as our basis for 
the metallicity gradient Z(R) and assume that both the elemental abundances and the gas 
phase abundances of all elements scale similarly with Z(R). Since both carbon and oxygen 
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are primary elements, the elemental carbon abundance gradient is seen to closely follow the 
oxygen gradient (Rolleston et al. 2000). We further assume that the grain size distribution 
does not vary significantly with Galactic radius, but that the total dust abundance scales with 
Z(R). Afherbach et al. (1997) note that they do not observe a jump in the nitrogen to oxygen 
abundance ratio at R < 6 kpc as suggested by Simpson et al. (1995). In addition, Rolleston et 
al. (2000) find that their data can be fit with a single slope of d[0/E]/dR = -0.067±0.008 dex 
kpc -1 without discontinuities (Twarog, Ashman, & Anthony-Twarog 1997) or a flattening of 
the gradient in the outer Galaxy as suggested by Fich & Silkey (1991). We take a constant 
gradient of d[0/H]/dR = —0.07 dex kpc -1 in the range 3 kpc < R < 18 kpc. This slope 
corresponds to a radial exponential scale length of = 6.2 kpc. 



3.4. Distribution of Ionized Gas 

Although the ionized gas in the Galaxy does not enter into our analysis directly, we 
include a brief discussion of it for completeness. Most of the ionized gas in the Galaxy is 
produced by photoionization. Taylor & Cordes (1993) identified three components of ionized 
gas: a diffuse component that extends out to R > 20 kpc, an annular component in the inner 
Galaxy centered at 3.7 kpc, and a component associated with spiral arms. Heiles, Reach, 
& Koo (1996) suggested that there are actually only two separate components, since the 
annular component is most likely due to spiral arms in the inner Galaxy. 

Taylor & Cordes (1993) infer that the diffuse component has a mean electron density 
in the Galactic plane at the solar circle of 0.019 cm" 3 , a vertical scale height of 0.88 kpc, 
and a radial distribution proportional to sech 2 (i?fc/20). Assuming that the ionization of He 
is similar to that of H (Slavin, McKee, & Hollenbach 2000), this corresponds to a surface 
density of diffuse ionized gas of Ennuis = 0.89sech 2 (i4/20) M pc~ 2 (not including He); 
at the solar circle, this is 0.75 M pc -2 . In the Wolfire et al. (1995b) model of the Galactic 
halo, the surface density of H in the collisionally ionized, hot interstellar medium (HIM) at 
the solar circle is 0.26 M pc~ 2 , which accounts for about 1/3 of the total diffuse ionized 
gas. We note that the HIM, however, has a much larger vertical scale height. 

The photoionized gas in spiral arms is associated with the H II regions produced 
by OB associations. These H II regions typically have dense cores, which appear as ra- 
dio H II regions, and much lower density envelopes, which absorb a significant fraction 
of the ionizing photons (Anantharamaiah 1985; McKee & Williams 1997). Most of the 
mass of the photoionized gas is in the envelopes, which have a surface density Shii,ciiv = 
3.5exp(-i4/3.5) M pc~ 2 (McKee & Williams 1997), or 0.31 M pc~ 2 at the solar circle. 
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4. Distribution of Energetic Photons and Particles 



The distributions of far-ultraviolet (FUV, 6 eV < hv < 13.6 eV) radiation, extreme 
ultraviolet (EUV, 13.5 eV < hv < 100 eV) radiation, soft X-ray (100 eV ^ hv ^ 1 keV) 
radiation, and cosmic rays are required to calculate the ionization fraction in the gas, the 
charge on grains, and the grain photoelectric heating rate. In the next three subsections we 
describe our adopted distributions. 



The FUV radiation field strength enters the photoelectric heating rate in two ways. 
First, it provides the total photon energy available for gas heating, and second, it governs 
the grain charge and thus the efficiency at which photon energy is converted into gas heating 
(Watson 1972). We calculate the far-ultraviolet radiation field by carrying out a simple 
radiative transfer calculation in an inhomogeneous medium. Note that we shall be scaling 
the numerical results to the measured value in the solar neighborhood and thus we need 
calculate only the variation from the local value. In addition, we are concerned with the radial 
variation of the mean intensity in the Galactic midplane and not the variation with height 
above the plane. And finally, we note that the vertical scale height of the H I and diffuse 
dust distribution, iff 1 , is always much greater than that of the OB stars that contribute to 
the FUV field. With these considerations, we assume that the FUV emissivity and opacity 
have a constant value between the midplane and the scale height of OB stars above the 
plane H° B (R), and that the emissivity is zero at heights greater than H° B (R). The mean 
intensity in the Galactic midplane at radius R is given by 



where b and / are Galactic latitude and longitude and J FUV (i?, I, b) is the far- ultraviolet 
intensity in direction (/, b) at Galactocentric radius R. The expression for the intensity in 
direction (/, b) is given by 



where j and k are the FUV emissivity and opacity respectively and r max is the maximum 
optical depth along the line of sight (l,b) (see the discussion below eq. [16]). The FUV 
emissivity is determined by the distribution of OB stars in the Galactic plane while the 
opacity is provided by dust, mainly in the diffuse atomic phases. Normalized to the local 
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value in the solar neighborhood, the FUV emissivity scales as 



j(Ro) 



j{R) 



SobCR) H? B (R ) 
Zob(Ro) H°*{R) 



(15) 



where Eob(-R) is the surface density of OB stars. We assume the vertical height H° B (R) is 
given by the H 2 height in the inner Galaxy i.e., H® B (R < R ) = 59 pc, and that it scales 
with the H 2 height at R greater than R . 

To proceed further we need to find how Eob(-R) scales with R, or equivalently, we require 
the radial exponential scale length, H^ B , of the OB star surface density. McKee & Williams 
(1997) examined the Smith, Biermann, & Mezger (1978) catalog of giant radio H II regions 
and found the surface density of the exciting OB associations between 3 kpc < R < 11 kpc 
can be fit with a scale length of roughly 3.5 kpc. They found no evidence for giant radio 
H II regions beyond 11 kpc and within 3 kpc (other than at the Galactic Center). Bronfman 
et al. (2000) recently obtained a scale length of only ~ 1.8 kpc for the surface density of 
embedded OB stars in the outer Galaxy (8.5 ~ R ~ 17 kpc). We note that the Bronfman 
et al. (2000) data can be fitted well by an _£^ B = 3.5 kpc scale length between 8.5 kpc and 
~ 13 kpc, with a steeper slope at greater radii. The Bronfman et al. results are based on 
observations of far-infrared radiation from deeply embedded OB stars, but the embedded 
stage may last a shorter time in the outer Galaxy because molecular clouds are smaller there: 
Solomon, et al. (1987) found no GMCs at R > 10 kpc (after correcting their distances to a 
Galactocentric distance of 8.5 kpc); Heyer et al. (2001) found no clouds more massive than 
10 5 M Q at R > 11.6 kpc, whereas half the molecular gas inside the solar circle is in clouds 
with M > 1O 6 M according to Williams & McKee (1997). Thus, it is conceivable that the 
total OB star distribution extends beyond 13 kpc without a break in slope. 

What do observations of other galaxies tell us about the radial distribution of OB stars? 
In M31, Cuillandre et al. (2001) have found B stars out to ~33 kpc, or 1.7 x R 25 , a distance 
corresponding to ~ 21 kpc in the Galaxy. Further evidence for star formation at large 
galactocentric distances is provided by Wang, Hoflich, & Wheeler (1997), who found that 
the radial distribution of SNe in disk galaxies is exponential (with an average radial scale 
length of 3.5 kpc) with no evidence for an outer cutoff in the distribution. Ferguson et 
al. (1998) found that star formation in several disk galaxies extended out to at least 2i? 25 
(corresponding to about 24.5 kpc in the Galaxy). They did find a break in the star formation 
rate, but the scale-height in the outer parts of the disk (i? 25 < R < 2_R 25 ) averages about 
0.3i? 25 , similar to the value found by McKee & Williams (1997) for the Galaxy between 3 
and 11 kpc. The H II regions in the outer parts of these galaxies were substantially smaller 
than in the inner regions, consistent with the lack of giant H II regions in the Galaxy beyond 
11 kpc. On the other hand, Martin & Kennicutt (2001) analyzed a larger sample of galaxies 
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and found that most galaxies exhibit a strong cutoff in their star formation at a radius 
determined by the Toomre criterion; this cutoff radius is generally of order i? 25 . In view 
of the uncertainties in the distribution of star formation at large distances, we shall extend 
our analysis only out to 18 kpc, corresponding to about 1.5i?25- Between 4 and 18 kpc, we 
shall adopt a scale length for the OB star surface density of H^ B = 3.5 kpc, noting that 
the actual distribution beyond 13 kpc is uncertain and, if anything, is less than our adopted 
distribution. Between 3 kpc and 4 kpc, we shall adopt a constant OB-star surface density; 
since the H I and H 2 surface densities appear to drop towards the Galactic center in this 
range, it is unlikely that the OB star surface density would continue to rise at radii less than 
4 kpc. 

Using Hft B = 3.5 kpc (at R > 4 kpc) and substituting iff 2 for the OB star vertical 
height, equation (15) can be simplified to 

j(Ro) Hf*(R) ' 1 j 

With the adopted cutoff in the OB star distribution, the maximum path of integration r max 
in equation (14) extends to R — 18 kpc in the plane or until the line of sight reaches the 
perpendicular height of z = H° B (R') = H° B (R )[H? 2 (R')/H? 2 (R )}. 

The opacity to the FUV radiation /tpuv(-R) depends on the dust abundance, which in 
turn is proportional to the mean gas density and metallicity. The run of mean density (n(R)) 
and metallicity Z(R) are discussed in §§ 3.1 and 3.3 respectively. Dust in the molecular 
phases is to a large extent shielded from the FUV radiation. Recent Orpheus I observations 
(Dixon, Hurwitz, & Bowyer 1998) confirm the earlier Copernicus result that the H to H2 
transition in diffuse gas occurs at E(B — V) ~ 0.1. For Ry = 3.1 this corresponds to 
Ay ~ 0.3 or Apuy ~ 0.6. Thus, the FUV field responsible for gas heating does not penetrate 
deeply into the molecular layer. In addition, the small volume filling factor of molecular gas 
(/ ~ 0.1%) means that the molecular component of the FUV opacity can be safely neglected. 

The opacity at visual wavelengths is found from the local observed extinction [Ay = 
N/(2 x 10 21 cm -2 )] and mean H I density ((nnif-Ro]) — 0.57 cm -3 ) from Dickey & Lockman 
(1990), yielding k v (R ) = 0.88 kpc -1 in the midplane. The FUV opacity for photoelectric 
heating is taken from the photodissociation region models of Tielens & Hollenbach (1985) 
and is based on the radiation transfer results of Flannery, Roberge, & Rybicki (1980). We 
use k fuv (R ) = 1.8k v (R ) = 1.6 kpc -1 . The FUV opacity is then given by 

KFW(R) = 1.6^5^> e -(^)/H- kpc -l 

with the radial metallicity scale length = 6.2 kpc. Since the scale height of diffuse gas is 
much greater than that of OB stars, we can consider the gas density in our integrations to 
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be independent of height and a function of only the radial distance R. Figure 4 shows the 
FUV opacity as a function of R. 

Locally the FUV (6 eV < hv < 13.6 eV) intensity has a measured strength of approx- 
imately 4irJ(R ) = 2.7 x 1(T 3 erg cm -2 s _1 (Draine 1978). This is a factor of 1.7 higher 
than the integrated field of Habing (1968), often used as a unit (G = 1) of flux in models 
of photodissociation regions. Since the local FUV flux is observed, we find the distribution 
of flux in the Galaxy by scaling j(Ro) in equation (14) so that the value of the integral at 
R = Rq is given by the Draine field. Results are shown in Figure 5. The drop off in intensity 
near 18 kpc is due to the abrupt cut off in the OB star population that we have imposed in 
the outer Galaxy. The true intensity distribution should fall off more gradually. 

We find that the calculated mean intensity at R > 4 kpc can be fitted by an exponential 
with a radial scale length of = 4.1 kpc, 

4ttJ fvy (R) = 47Tj FVV (R )e-^ R ^ H 'n J R>4kpc, (18) 

whereas for R between 3 kpc and 4 kpc it is constant, 

4ttJ fuv ( j R) = 47rJ FUV (i? )e- (4 - i?o)/H « 3 > R < 4 kpc. (19) 

These expression are good to within 10% for R between 3 kpc and 17 kpc. 

4.2. Cosmic Rays 

The low-energy (E < 100 MeV) cosmic rays that contribute to ionizing the CNM and 
WNM do not travel far from their point of origin (Kulsrud & Cesarsky 1971; Spitzer & 
Jenkins 1975). We obtain the cosmic ray ionization rate as a function of position in the 
Galaxy by scaling the local primary rate (taken to be (qr = 1-8 x 10~ 17 s^ 1 ) by the production 
rate of cosmic rays per unit area (sources) divided by the mass per unit area (sinks). For the 
distribution of sources we use the surface density of OB stars (§ 4.1), and for the distribution 
of sinks we use the total surface density of molecular and neutral atomic (CNM+WNM) gas 
(§§ 3.1 and 3.2). This scaling differs from that adopted by Hunter et al. (1997) and Bertsch 
et al. (1993), who made the assumption that the cosmic ray intensity is proportional to 
the surface density of gas alone; however, these authors were studying Galactic gamma-ray 
emission, and were therefore interested in cosmic rays with higher energies than those that 
dominate the ionization. Such high energy cosmic rays can travel more freely in the Galaxy 
and therefore acquire a fairly homogeneous distribution. 

The surface density of WIM can potentially affect our cosmic-ray ionization rate by 
providing an additional sink for cosmic rays (the energy loss rate for cosmic rays in ionized 
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gas is several times that in neutral gas — Ginzburg & Syrovatskii 1964). In order to simplify 
our model, we have chosen to neglect the effect of the WIM surface density in calculating the 
distribution of sinks of cosmic rays; note that this effect enters only insofar as the distribution 
of the WIM differs from that of the rest of the gas. Our neglect of the effect of the WIM 
will not strongly influence our results since EUV/X-ray radiation generally dominates the 
ionization and FUV radiation dominates the heating. 

The resulting primary cosmic ray ionization rate is given by: 
Ccr = 1.8 x i -ir e -(R-Ho)/HOB 

where H° B = 3.5 kpc, S ncutral (i?) = S H2 (i?) + S H i(i?) and S neutral (i? ) = 6.4 M pc" 2 . The 
variation in the primary cosmic ray ionization rate with Galactocentric radius R is shown in 
Figure 6. 

4.3. Soft X-ray and EUV Radiation 

As demonstrated in WHMTB, the ionization of H and He by soft X-rays and EUV 
photons (13.6 eV< hv ^ 10 3 eV) provides a source of electrons in the local WNM amounting 
to an electron fraction of n e /n ~ 2% (at n = 0.3 cm" 3 ), a fraction that depends on the column 
N w of warm absorbing gas and dust traversed by the X-rays. These electrons neutralize 
the positive charging of grains caused by the FUV photoelectric effect, and thereby help 
maintain a high photoelectric heating efficiency in the WNM phase. The local X-ray field 
arises from the Local Bubble and the Galactic halo plus an extragalactic background, with 
the low energy (50 eV < E< 100 eV) emission from the Local Bubble dominating the 
ionization at typical columns N w pa 10 19 cm" 2 . WHMTB used a fit to the observed X-ray 
intensity from Garmire et al. (1992) to generate the X-ray spectrum incident on a WNM 
"cloud" . The fit consisted of the temperature and emission measure in the Local Bubble and 
halo components, an absorbing layer through which the halo is seen, and an extragalactic 
component and absorbing layer. The ionization rate at a column N w into the WNM cloud 
was then calculated by computing the attenuation of the incident X-ray flux by the column 
N w . 

The Garmire et al. (1992) fit was based on pre-ROSAT data. Snowden et al. (1998) and 
Kuntz & Snowden (2000) evaluated the ROSAT data to analyze the origin and distribution 
of the diffuse X-ray background. They find somewhat lower temperatures than Garmire et 
al. (1992) for the Local Bubble and argue that the 1/4 keV and 3/4 keV band observations 
can be fit only if the halo is emitting at two or more temperatures. Using the observed ratios 
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in the ROSAT Rl, R2, R3, and R4 bands, Kuntz & Snowden (2000) find log T b = 6.11 for 
the bubble and log = 6.06, log T^ 2 = 6.46 for the two halo temperatures as compared 
with log TJ, = 6.16 and log = 6.33 from Garmire et al. (1992). (Note that there is an error 
in Table 1 of WHMTB in that the emission measures for the bubble and halo components 
are reversed). 

The fit to the ROSAT data provides the local soft X-ray flux. However, we know that 
the Sun is located in a bubble of hot gas that dominates the low energy radiation and hence 
the ionization, and the local emission may differ from the average emission at the solar circle. 
We derive the average X-ray emission in the disk from the calculations of Slavin, McKee, & 
Hollenbach (2000). They define the emissivity per unit area {e V A) as 



where (e u ) is the volume emissivity and the integration is carried out perpendicular to the 
disk. The average value of the mean intensity is found from (e v A) as 



where r v ^ is the opacity through the disk, and r\ v is the mean escape probability of soft 
X-rays out the Galaxy. The opacity r v ^ is taken to be that of a cloudy medium of column 
density through the disk A^ WNM ,d with clouds of typical column density iV cl (H I). (These 
parameters appear as t 0u , Nho±, and iVno c respectively in the Slavin et al. notation). We 
consider WNM clouds only; the small filling factor of the CNM means that it does not 
contribute much to the opacity of very soft X-rays. The opacity is then given by t v ^ = 
[iVwNM,d/-Wci(H I)][l — exp(— 7V )C i)] where t VjC \ is the optical depth through a single WNM 
cloud, and iV WN M,d is the WNM column through the disk. 

Slavin et al. determine the average emissivity produced by supernova remnants (e^ NR ), 
which accounts for the time and space averaged emissivity produced as supernova remnants 
evolve into the interstellar medium of ambient density n a , although Slavin et al. show that 
the average emissivity is very insensitive to n a . Slavin et al. applied their model to a specific 
line of sight towards a high latitude cloud and calculated the ionizing flux as a function 
of height in the disk and found that they could successfully match the observed fractional 
ionization. 

In addition to the soft X-rays emitted by the supernova remnants, stellar EUV radia- 
tion also contributes to the ionization rate in the WNM. We use the stellar EUV spectrum 
shown by Slavin & Frisch (1998), which is derived from Extreme Ultraviolet Explorer ob- 
servations (Vallerga 1996) and corrected for extinction by the local interstellar cloud (taken 
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to be iV(HI) = 9 x 10 17 cm' 2 .) While keeping the stellar EUV spectral shaped fixed, we 
adjust the level of the stellar EUV emissivity (e*) so that the total emissivity per unit area 
( € a) = / du dz [(e^ NR ) + (e*)] provides the observed ionizing photon flux of the Galactic disk 
outside H II regions, as deduced from Ha observations ((e^) ~ 4 x 10 6 phcm' 2 s -1 ; Reynolds 
1984, 1995). We find approximately 43% of the EUV photon flux comes from stars; this is 
consistent with the results of Slavin et al. (2000), who suggest that hot gas in SNRs produces 
about half the total number of ionizing photons in the diffuse ISM. 

We adopt parameters for the total column density of WNM through the disk appropriate 
for the Galactic average at the solar circle. Dickey & Lockman (1990) fit the vertical distri- 
bution of H I in the solar neighborhood with three distinct components. Following Kulkarni 
& Heiles (1987), we identify the two components with the largest scale heights with the 
WNM. The total WNM column density through the disk is then N WNM ,d(Ro) = 3.45 x 10 20 
cm' 2 , or ~ 56% of the total H I column at R . We also set the mean escape probability 
equal to zero {r\ v = 0). For our standard model we set iV cl (H I) = 1 x 10 19 cm' 2 , comparable 
to the column densities of the WNM clouds observed along the line of sight toward the halo 
star HD 93521 (Spitzer & Fitzpatrick 1993) and disk star 7 2 Vel (Fitzpatrick & Spitzer 1994; 
this line of sight contains 4 WNM clouds ranging in H I column density from 2.7 x 10 18 cm' 2 
to 4.5 x 10 19 cm -2 with an average value of 1.5 x 10 19 cm -2 ). Note that in our formalism, N c \ 
is the WNM column density in a typical cloud that provides the opacity for the EUV and 
soft X-ray radiation. Note also that the phase diagrams and thermal processes presented in 
§ 6 apply to the WNM/CNM boundary within a cloud and thus radiation incident upon the 
cloud must pass through an additional column N d (or N w in WHMTB notation). 

Using the flux at z — 0, we find that at a cloud column of AT cl (H I) = 10 19 cm' 2 , the 
primary EUV plus X-ray ionization rate of hydrogen is Cxr = 1-6 x 10' 17 s' 1 , a factor 
~ 1.6 lower than that used by WHMTB. (This rate is approximately equal to the primary 
cosmic ray ionization rate, 1.8 x 10' 17 s' 1 . The total ionization rate from either cosmic 
rays or EUV/soft X-rays is larger than the primary rate due to the effects of secondary 
ionizations. The secondary rate increases with the energy of the primary ejected electron 
and with decreasing ionization fraction. In the WNM at a cloud column of iV c i = 10 19 cm' 2 
and density n ~ 0.3 cm' 3 the total EUV/X-ray rate is ~ 5.3 x 10' 17 s' 1 , about 1.5 times 
higher than that from cosmic rays, while in the CNM at a density of n ~ 33 cm -3 the total 
EUV/X-ray rate is ~ 7.5 x 10' 17 s' 1 , or ~ 2.7 times higher than that from cosmic rays.) 
The ionizing photon intensity incident on clouds is ~ 9.4 x 10 3 cm' 2 s' 1 sr' 1 . This is a 
factor 2.2 larger than that obtained by Slavin et al. (2000) because (1) we are modeling a 
typical region of the ISM, which has a significant flux of stellar ionizing photons and (2) we 
have adopted a somewhat smaller cloud column density. 
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Having determined the value for the soft X-ray ionization rate at R = Rq, we must 
now scale the results for other Galactic radii. The ultimate energy source for the hot gas 
that produces the X-ray flux is supernova explosions, while the opacity arises from the 
surface density of (WNM) H I gas. Thus, we assume that the X-ray ionization rate per 
hydrogen in the Galactic midplane scales as the OB star distribution divided by the H I 
surface density, Cxr(-R) oc £ob(-R)/£hi(-R), as shown in Figure 6. We assume the stellar 
EUV photoionization rate scales the same as the soft X-ray ionization rate. Because of the 
increased gas opacity to the EUV/X-ray radiation compared to FUV radiation, the mean 
free path for the EUV/X-ray photons are much shorter than for FUV. Thus, the dust in 
CNM clouds dominates the FUV opacity, but the gas in WNM clouds dominate the EUV/X- 
ray opacity. Since the fraction of WNM column to the total is unknown outside the solar 
neighborhood we have assumed, for the purposes of estimating the X-ray distribution, that 
the fraction is always given by the value in the solar neighborhood. If the WNM fraction were 
to increase in the outer Galaxy, for example, the X-ray ionization rate would be lower than 
our adopted rates. We also include ionization by the extragalactic X-ray and EUV radiation 
field of Sternberg, McKee, & Wolfire (2002) which is based on the work of Haardt & Madau 
(1996) and Chen, Fabian, & Gendreau (1997). The extragalactic field passes through an 
absorbing column given by one half of the total WNM column density. At the solar circle 
the extragalactic field provides ~ 1% of the total ionization rate (for N c \ = 10 19 cm" 2 ), while 
at it! = 18 kpc the extragalactic rate rises to ~ 12% of the total rate. 



5. Chemical and Thermal Processes in the Neutral Phases 

The chemical and thermal processes included in this work are slightly modified from 
those discussed in WHMTB. The main changes involve the PAH reaction network and PAH 
reaction rates. These rates are important because reactions with PAHs affect the PAH charge 
state and electron abundance which in turn affects the photoelectric heating rate. We have 
simplified the PAH chemical reaction network compared to that used in WHMTB. Here we 
use rates appropriate for a single PAH size containing Nq = 35 carbon atoms, the mean size 
in the distribution between 3 A and 15 A. We have also dropped the adsorption reactions 
used in WHMTB since these were found to be not important in determining the PAH charge. 

For the PAH reaction rates we use the photoionization/photodetachment rates of Bakes 
& Tielens (1994) and a modified form of the Draine & Sutin (1987) formalism for the inter- 
action between ions and electrons with neutral and charged PAHs in which we multiply all 
of the collisional rates by a factor 0pah = 0.5. There is considerable uncertainty in applying 
a classical treatment of the interaction between atoms and grains to the molecular regime. 
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In particular, collision rates (using 0pah = 1) consistently overestimate the laboratory mea- 
sured rates for electron attachment to neutral PAHs and for electron recombination with 
PAH + (Allamandola, Tielens, & Barker 1989; Salama et al. 1996; Weingartner & Draine 
2001a and references therein). We will consider 0pah to be a parameter and rely on obser- 
vation to guide us in its appropriate value. As noted by Lepp et al. (1988) and recently by 
Bakes & Tielens (1998), Welty & Hobbs (2001), and Weingartner & Draine (2001b), ion re- 
combination on small grains and PAHs can be important in determining the neutral fraction 
of metals. (In Appendix C we also elucidate the conditions under which reactions of C + and 
H + with PAHs can affect the electron abundance.) Specifically, we find that the abundance 
of neutral carbon is sensitive to the rate of C + recombination with PAH - . We find that the 
observed C I/C II column density ratio in diffuse clouds ( ^ 3x 10~ 3 ; Welty & Hobbs 2001; 
Jenkins & Tripp 2001) imply 0pah — 0.5 and we adopt this value for our standard model. 
In § 6.1 we discuss the effect of higher and lower values of 0pah on the phase equilibrium of 



We have also updated the gas phase chemical reaction rates according to the list of 
Millar, Farquhar, & Willacy (1997) and modified the H2 formation and dissociation rate 
according to the discussion in Kaufman et al. (1999). These additional chemical changes 
mainly affect the molecular pathways included in our network and are of minor consequence 
for the atomic phases discussed in this paper. Note, however, that the adopted gas phase 
carbon and oxygen abundances have been modified from their WHMTB values (§ 3.3). The 
results of these changes are discussed in §§ 6.1 and 6.2. 

The dominant heating process at R = Rq in both the CNM and WNM phases is grain 
photoelectric heating. We use the heating rate determined by Bakes & Tielens (1994), 
modified by the higher PAH abundances as explained in § 3.3, and modified by the parameter 
0pah which scales the electron- PAH collision rates. For the same density, electron abundance, 
temperature, and incident FUV field, the higher abundances result in an increase in the 
heating rate by a factor 1.3. The cooling rate due to electron recombination with PAHs 
increases similarly by a factor 1.3. In Bakes & Tielens (1994) the heating/cooling rates were 
calculated self consistently with the PAH ionization state as a function of the photo and 
collision rates. To maintain this self consistency, the heating/cooling rates (and the fit to 
the rates) must be modified for the parameter 0pah- Including the correction for higher PAH 
abundances, the heating rate per unit volume is given by 
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where n is the hydrogen nucleus density, and the heating efficiency e is given by 

4.9 x 10- 2 3.7 x 10- 2 (T/10 4 ) - 7 
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The cooling rate per unit volume is given by 



n 2 A = 4.65 x Kr 30 T' 




n e 4>PAKn erg cm 3 s 



i 
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with (3 = 0.74/T - 068 . 

The dominant cooling process in the CNM is radiative line cooling in the [C II] 158 jum 
fine-structure transition. Cooling in the WNM is provided by several processes: radiative 
line cooling by [C II] 158 fim, [O I] 63 jum, and Lya, and by electrons recombining onto 
grains (refer to details in WHMTB.) In this paper we adjusted the collision rate for the 
excitation of C + by impacts with e~ to that of Blum & Pradhan (1992). 



Utilizing the results from the previous sections on the FUV intensity and ionization 
rates in the Galaxy as functions of R, we calculate phase diagrams - gas thermal pressure 
P versus hydrogen nucleus density n - for several Galactocentric distances R. The curves 
are generated by calculating, at constant n, the chemical equilibrium abundances and the 
thermal equilibrium temperature. We then step through n and plot the calculated thermal 
pressure, P = Sn^/cT, where % ranges over all chemical species. The calculations are carried 
out for various WNM cloud columns N d and apply to the WNM/CNM boundary. An 
appropriate range is 3 x 10 18 cm -2 < N c \ < 1 x 10 20 cm -2 . The upper limit is set by the 
size scale at which turbulent pressure begins to dominate. As discussed in § 2, we obtain 
£p(WNM)~ 215 pc, or 7V cl ~ 2 x 10 20 cm~ 2 for a typical density n ~ 0.3 cm~ 3 . The lower 
limit is set by the requirement that there be a substantial neutral fraction in the cloud. We 
note that the smallest neutral column density in the warm clouds along the line of sight to 
the disk star 7 2 Vel is 2.7 x 10 18 cm" 2 (Fitzpatrick & Spitzer 1994). The column density of 
the ionized clouds in the MO model of the ISM is 2.2 x 10 18 cm -2 ; since this model agrees 
with observations of Ha and pulsar dispersion measures in the disk, it is difficult to have 
predominantly neutral clouds that are smaller than this. We shall use N c \ = 10 19 cm -2 as 
a standard and demonstrate the effects of both higher (1 x 10 20 cm~ 2 ) and lower (3 x 10 18 
cm~ 2 ) columns. 

We present in Figure 7 phase diagrams for Galactic radii R = 3, 5, 8.5, 11, 15, and 
18 kpc. The model parameters are given in Tables 1 and 2. As demonstrated by Field 
(1965), the region of thermal stability in P versus n phase diagrams lies in the range where 



6. 



Results 



6.1. Phase Diagrams 
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dP/dn > 0. Where dP/dn < the gas is thermally unstable to isobaric perturbations. If 
the pressure curve has a characteristic shape shown in Figure 7, two thermally stable phases 
may coexist in pressure equilibrium within a range of gas pressures, P min to P max . At thermal 
pressures greater than P max only the cold phase (CNM) may be present, while at thermal 
pressures less than P min only the warm phase (WNM) may exist (see also further discussions 
in reviews by Shull 1987 and Begelman 1990). 

We see that a two-phase (WNM+CNM) equilibrium is possible in the Galactic midplane 
at all Galactic radii between 3 kpc and 18 kpc. The pressure range P m i„ and P max are listed 
in Table 3 for each radius along with the range in gas temperature and density for WNM 
and CNM gas between P min and P max . Also listed is the average pressure, where we adopt 
the geometric mean of P min and P max to represent the average, 

Pth,ave(P) = [Pmin(P)Pmax(P)] 1/2 • (26) 

We also give the temperature and density at Pth, av e for WNM and CNM gas. Compared to 
the results in WHMTB at R = 8.5 kpc, P m i n is higher by a factor of ~ 2.0 and P max is higher 
by a factor of ~ 1.3. The difference in P m i n is mainly due to the revised (lower) gas phase 
carbon and oxygen abundances from those in WHMTB. The lower abundance of coolants 
results in a higher gas temperature. The higher P max is partly due to the lower abundance 
of coolants, but mitigated by the effects of the collision rate parameter 0pah- The collision 
parameter affects the electron fraction, the photoelectric heating rate, and the cooling rate 
due to electrons recombining onto positively charged grains. Lower values of 0pah result in 
a higher electron fraction in WNM gas. At P max , higher electron fractions and enhanced 
recombination cooling plays a role in limiting the maximum pressure. 

We next examine the effects of the PAH collision rate parameter 0pah and the PAH 
abundance on the pressure curves. We show in Figure 8 and list in Table 4 results for R = 8.5 
kpc and 0pah = 0.25, 0.5 and 1.0, with the standard PAH abundance (^pah / n = 6 x 10~ 7 or 
an amount of C in PAHs of 22 x 10 -6 relative to H). We also show results for 0pah = 0.5 and 
a lower PAH abundance of npAu/n = 4x 10~ 7 as used by WHMTB. Results for Go = 1-1 as is 
appropriate for the interstellar radiation field of Mathis et al. (1983) are also given in Table 4 
but not shown in the figure since the resulting pressures are very similar to the case for low 
PAH abundance. Over the range of 0pah from 0.25 to 1, P m i n changes by a factor of only 1.5 
while P max changes by 1.9 mainly due to the increased effects of recombination cooling. We 
conclude that the results for the phases of the ISM are very robust against variations in the 
PAH physical and chemical characteristics. For example, the average temperature hardly 
varies. The largest change occurs for 0pah = 0.25 in which P max decreases by ~ 35%. 

The phase diagrams presented in Figure 7 used a column of atomic gas N c \ = 10 19 cm -2 . 
We illustrate in Figure 9 the effects of higher (N d = 10 20 cm" 2 ) and lower (N d = 3 x 10 18 
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cm~ 2 ) column densities at R = R = 8.5 kpc. As the column density decreases at a given 
Galactic radius, the electron abundance rises (due to hydrogen photoionization by EUV/X- 
rays). The increased electron abundance neutralizes the grains and enhances the rate of 
grain photoelectric heating. In addition, the EUV/X-ray heating rate rises. Therefore, lower 
column densities result in higher temperatures and pressures, P m i n and P max - At columns 
of N c i = 3 x 10 18 cm' 2 , P min /k = 2560 K cm" 3 and P max /A; = 7830 K cm" 3 while at the 
higher column of N cl = 1 x 10 20 , P min /k = 1240 K cm" 3 and P max /k = 2310 K cm" 3 . The 
calculated range in two phase (CNM+WNM) thermal pressures over the column densities 
N c \ = 3 x 10 18 cm' 2 to N c \ = 1 x 10 20 cm" 2 are in reasonably good agreement with the 
thermal pressures of P/k — 10 3 — 10 4 K cm -3 that are observed in the local ISM (Jenkins et 
al. 1983; Jenkins & Tripp 2001). For the fiducial case N c i = 10 19 cm -2 , the average pressure 
at the solar circle is predicted to be 3070 K cm -3 , consistent with the C I population ratio 
observed by Jenkins & Tripp (2001) (see § 8.1). 



6.2. Thermal Processes 

Figures 10a through d show the dominant thermal processes at Galactic radii of R — 5, 
8.5, 11, and 17 kpc. We see that at all radii, photoelectric heating dominates in both the 
CNM and WNM phases. The [C II] 158 jum line emission dominates the cooling in the CNM 
and is generally a factor of ~ 20 weaker (per H) in the WNM. The photoelectric heating 
rate per hydrogen is nearly the same as the results presented in WHMTB with the current 
rates approximately 20% lower at a density of 30 cm" 3 . The higher PAH abundance used in 
this paper is offset by a lower gas phase carbon abundance (and lower electron abundance in 
CNM gas), higher temperatures, and a 0pah = 0.5 collision rate parameter, which decreases 
the photoelectric heating efficiency Due to the higher CNM temperatures however, we arrive 
at an average CNM density which is a factor ~ 2 lower than WHMTB (see Table 4), which 
results in an average C II cooling rate per hydrogen a factor of ~ 1.5 lower than our previous 
result. 

It is instructive to consider the total cooling luminosity of the ISM, which is dominated 
by the [C II] 158 fim line. The [C II] luminosity of the Milky Way has been measured by 
COBE to be 5 x 10 7 L (Wright et al. 1991). Tielens (1995) considered the possible global 
sources of this emission. In the inner galaxy, much of it may originate in extended low 
density H II regions (Heiles 1994) while most of this emission in the Galaxy as a whole has 
to stem from the CNM (WHMTB). The WNM can not contribute much to the total [C II] 
luminosity due to the low emission rates in such tenuous gas. We consider the potential 
for mechanical heating by supernovae in Appendix B and conclude that neither the shock 
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heating nor the turbulence generated by SNRs can contribute substantially to heating the 
H I gas. (Of course, mechanical energy dominates the energy balance of the hot ionized 
medium and coronal gas in the Milky Way.) Thus, the thermal structure of the ISM is 
largely dominated through the coupling of the gas to the stellar, non-ionizing radiation field. 
The fraction of [C II] emission in the inner Galaxy which arises from CNM or warm ionized 
gas depends on the mass and density (or filling fractions) of these components. Fitting the 
observed profile of [C II] emission versus Galactic longitude may provide a test of the models 
presented here, however a detailed model of the ionized gas (in pressure equilibrium and in 
overdense regions) is beyond the scope of this paper. We shall discuss in a future paper the 
implications of our [C II] emission rates for the filling fractions of the WNM and CNM gas, 
and for the origin of the [C II] emission in the Galaxy. 



6.3. Predicted Infrared Radiation Field 

In this subsection we discuss a check on the distribution and local values of the model 
opacity and interstellar radiation field by comparing the calculated infrared intensity emitted 
by dust with observations from the COBE satellite. 

We calculate the integrated infrared continuum intensity along a line of sight assuming 
that the predominantly far-infrared emission is optically thin, 

I m = f dvl u = [ ds' f dv I ' daB v [T{a)\n v {a) , (27) 
Jm, Jo JiR J 

where B u [T(a)] is the Planck function for the grain-size dependent temperature T(a) and 

K v (a) is the grain absorptive opacity for grains of size a. In thermal equilibrium, this emission 

just balances the heating by UV photons, 



J IR = f ds' f dv f da Jl SRF 
Jo Juv J 



K v {a) , (28) 



where J* SRF is the interstellar radiation field. We have ignored the small amount of energy 
that goes into photoelectric heating of the gas. We take the spectral energy distribution of 
the interstellar field to be that of Draine (1978) between 912 A< A < 2000 A, van Dishoeck 
& Black (1982) between 2000 A< A < 3400 A, and Mathis et al. (1983) for A > 3400 A. 
The grain absorptive opacity k v is calculated as in Wolfire & Cassinelli (1986) and Wolfire & 
Churchwell (1994) using a grain abundance and optical constants from Draine & Lee (1984, 
1987) with grains distributed in size as a power law (n[a\ da oc a~ 3 - 5 da; Mathis et al. 1977). 
The dust opacity is scaled with R in proportion to the mean H I density and metallicity, 

k v {R) = k u (R ) j nH j^. e-^-^ (29) 
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while the interstellar radiation field is scaled with R according to results reported in § 4.1. 

We show in Figure 11 the calculated IR emission compared with the COBE longitudinal 
profile reported by Sodroski et al. (1994). We obtain a reasonably good fit. We have not 
compared with the region within 3 kpc of the Galactic Center, since we have not modeled 
the H I there. The largest deviations occur at the locations of the spiral arms, where some of 
the FUV radiation is absorbed by molecular gas associated with star-forming regions; recall 
that we have not included this in our models, which are for the diffuse ISM. Integrating over 
galactic longitude, but excluding the region near the Galactic Center, the mean intensity 
from COBE is 1.4 x 10~ 5 W m -2 sr^ 1 , whereas our model gives 1.1 x 10~ 5 W m -2 sr _1 . 
According to Parravano et al. (2002), about 80% of the FUV radiation is absorbed in the 
diffuse ISM, and 20% is absorbed by the nearby natal giant molecular cloud. Using the 
interstellar radiation field and opacity given in the preceding paragraph we estimate that 
approximately half the dust heating is produced by the FUV radiation from OB stars and 
half by longer wavelength radiation from older stars generally located far from molecular 
clouds. Correcting for the OB starlight absorbed by nearby molecular clouds, our model 
gives 1.2 x 10~ 5 W m -2 sr" 1 for the total IR emission, in reasonably good agreement with 
observation. 



We have presented phase diagrams which show P min and P max , and have presented 
observational evidence that both CNM and WNM gas exists, but we have not yet tried to 
estimate the thermal pressure in order to make comparison to P m i n or P max - One way to 
determine the pressure is from the density and volume filling factors of the H I gas. Our 
results are in terms of the local density n, but they can be recast in terms of the mean atomic 
hydrogen density (uhi) which is often better determined observationally. Let nwNM(-Pmin) be 
the equilibrium density of the WNM at the pressure P mm , etc. We generalize the treatment of 
Krolik, McKee, & Tarter (1981) (who used the notation n for the mean density) by allowing 
for the possibility that the two-phase medium fills only a fraction /hi of the volume and by 
allowing for pressure fluctuations in the gas. We follow Krolik et al. in assuming that the 
gas is in thermal equilibrium, which is not strictly true in a turbulent medium. 

The ratio of the CNM mass to the WNM mass is (Krolik, McKee, & Tarter 1981) 



6.4. Multiphase ISM 



MqNM 
M W NM 



[(»hi)/(»wnm/hi)] - 1 
1 - [(^hi)/(^cnm/hi)] ' 
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Similarly, one can show that the ratio of the volume filling factors is 



/cNM 
/wNM 



1 - (n WNM /Hi/ (^hi)) 
("cnm/hi/ ("hi)) - 1 ' 



(31) 



Let 1Z = ncNM(-Pmin)/ n WNM(-Pmax); if 7Z is large, then there is an extensive range of density 
in which the H I is thermally unstable. The cases summarized in Table 3 all have TZ > 6. 

There are five regimes for two-phase media: 

1- (ft-Hi) //hi < ^WNM(-fmin)- Since the average density of the H I is less than the minimum 
required for CNM to exist, it follows that most, if not all, of the volume of H I is filled 
with WNM. If the gas is isobaric, the H I is all WNM. 

2. n WN M(-P m in) < (^hi)//hi < %nm (-Pmax) • Here again, most of the volume of H I must 
be filled with WNM (one can show from eq. [31] above that /cnm//wnm < !])> 
but even for the isobaric case some CNM can coexist with the WNM. In the presence of 
turbulent pressure fluctuations, some WNM will be driven into CNM, and both phases 
will occur 

3. nwNM(-Pmax) < (nm) / fni < ^cnm ( -Pmin ) • If the local density n were in this regime, 
the gas would be thermally unstable and two phases would form. As a result, in this 
case two phases must exist. 

4- ^CNM(-Pmin) < (^hi)//hi < (-Pmax) • 111 this regime, most of the mass must be 

CNM (one can show from eq. [30] above that Mcnm/^wnm > TZ — 1), but even for the 
isobaric case some WNM can occur. In the presence of turbulent pressure fluctuations 
(which are negative as well as positive), some CNM will be driven into WNM, and 
both phases will occur. However, in this regime, it is possible that all the neutral gas 



5. ^CNM(-Pmax) < (^hi)//hi- Finally, in this case most, if not all, the H I must be in the 
form of CNM. If the gas is isobaric, all the H I must be CNM. 

According to Table 3, the local ISM is characterized by nwNM(-Pmm) = 0.21 cm~ 3 , 

flWNM(fmax) = 0.86 CHI" 3 , 7l C NM(-Pmin) = 6 -9 CHI" 3 , and 7l C NM(-Pmax) = 71 CHI" 3 . The locally 

observed value of the mean atomic hydrogen density (nm) m the Galactic plane is 0.57 
cm -3 . Therefore, if the H I is pervasive, the local ISM would be in the regime in which 
most of the volume of H I is primarily WNM. Observations show that most of the mass of 
the H I in the local plane is CNM (Dickey & Lockman 1990). Recently, Heiles & Troland 
(2002), utilizing H I 21 cm measurements, estimated a local value of /hi = 0.5 based on 



is CNM. 
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the assumption that the local thermal pressure of WNM is 2240 K cm -3 and T WNM = 4000 
K so that «wnm = 0.56 cm~ 3 . In this paper, we argue that the local WNM pressure is 
approximately 3100 K cm -3 and that T WNM ~ 8000 K so that «wnm — 0.35 cm -3 . The 
Heiles & Troland result then becomes /hi — 0.79. These results lie close to the theoretical 
estimates which range from /hi = 0.4 (MO) to 0.8 (Slavin & Cox 1993). If a substantial 
fraction of the volume of the ISM is hot then /hi ~ 0.5 and the local ISM is marginally in 
the regime in which there must be two phases. 

We extend our analysis to the inner and outer Galaxy by estimating the thermal pressure 
in the midplane in the limiting case in which the H I is entirely in the form of WNM. We derive 
a WNM pressure and check to see if it is self consistent (i.e., is the thermal pressure less than 
P min or P max ) to assume that all the neutral gas is WNM. To calculate the thermal pressure 
we rely on the assumption that the total pressure in the midplane is balanced by the weight 
of the overlying gas layers (i.e., hydrostatic equilibrium). We first consider the case in which 
the Galactic H I resides in a continuous, thermally supported WNM layer, and then modify 
this thermal pressure calculation for the effects of non-thermal support, including cosmic 
rays, magnetic fields, and turbulence. We also estimate the thermal pressure by using the 
mean atomic hydrogen density (nni) shown in Figure 3 to derive a "mean" thermal pressure 
which we compare to P min and P max to find the region over which a two-phase medium must 
exist. 



6.4.I. Thermally Supported WNM 

For an isothermal gas in hydrostatic equilibrium in the Galactic gravitational field, the 
vertical density distribution is given by 

n(R, z) = n (P)e [ * ( *' 0) -* (H ' 2)1/< ^ (32) 

where n(R, z) is the gas density at Galactic radius R and height above the plane z, n (R) = 
n(R,z = 0) is the density in the midplane, z) is the Galactic gravitational potential, 
and <7 t h is the isothermal sound speed. We take <3?(-R, z) from a variant of model 2b of Dehnen 
& Binney (1998), which has a disk mass of 5.3 x 10 10 M (by comparison, the disk model 
used by Wolfire et al. 1995b had a mass of 1.0 x 10 11 M ). The numerical code to calculate 
the potential was kindly provided by W. Dehnen. Each density component in the disk is 
assumed to be of the form 

Yl d / R m R \z\\ , 

P=77-exp ' 33 

2z d \ R Rd z d J 

where z d is the vertical scale height and R d is the radial scale length of the disk. The 
parameter R m allows for the depression in the gas density observed in the inner several 
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kiloparsecs of the Galaxy; for the stellar disks, they set R m = 0. Dehnen & Binney (1998) 
assumed that the gas in the disk could be described by a single component with = 40 pc 
and, for Model 2b, = 5.1 kpc. A limitation of their model (which is relatively unimportant 
for their application) is that the vertical scale height for each component is assumed to be 
independent of radius. We have altered their model for the gas to make it more consistent 
with the discussion in § 3, with one component for the molecular gas, two components for 
the H I in order to capture the radial variation, and one component for the H II. (We have 
included only the diffuse H II component discussed in § 3.4, since the ionized gas contributed 
by the envelopes of H II regions is always very small compared to the column density of stars 
plus H I. We have approximated sech 2 [i? fc /20] as exp[— R k /30\). We have also altered their 
model for the thin stellar disk to allow the scale height to increase in the outer Galaxy. For 
Galactic radii R > Rq we assume that Zd oc 1/S where E includes the total (gas plus stars) 
surface density and then recalculate the potential for each radius R using the appropriate 
Zd- The parameters entering the fit are summarized in Table 5. The total surface density 
at the solar circle is 10.1 M Q pc~ 2 . This is 2.5 M Q pc~ 2 less than that of Dehnen & Binney 
(1998), so we added 2.5 M Q pc~ 2 to the thin stellar disk in order to maintain the same value 
of the total surface density. 

We wish to calculate the pressure that the WNM in the Galactic plane would have 
if all the H I in the disk were in the form of WNM and if the pressure support of the 
WNM were entirely thermal. We label this hypothetical pressure Pwnm'- It is given in 
terms of the midplane H I density ra H i,o(-R) by Pwnm'(R) = l-lnm jQ (R)kT. To determine 
n ni,o(R), we first calculate the column density from equation (32), noting that the the 
velocity dispersion of the neutral WNM is <j t h = 7.2 km s _1 (T = 8000 K). We then equate 
this theoretical column density to the observed column density, Nm(R) = 2 J njn(R,z) dz, 
and solve for nni,o(-R)- Figures 12a through 12c show the pressure obtained in this manner 
for N c \ = 1 x 10 19 cm -2 , 1 x 10 20 cm -2 , and 3 x 10 18 cm~ 2 , respectively, compared to P m i n , 
P max , and the average thermal pressure Pth.ave(-R)- The kinks in Pwnm' are a reflection of the 
kinks in our adopted H I surface density (eqs. [7] through [10]). In Appendix A we provide 
an analytic solution for Pwnm' at the solar circle which provides a value of Pwnm' = 7800 
K cm -3 , in good agreement with the numerical solution of Pwnm' = 8615 K cur 3 . These 
values exceed P max for iV c i ^ 3.0 x 10 18 cm -2 (see Table 3), so we conclude that a thermally 
supported H I layer cannot be all in the form of WNM gas and some of the H I must be 
forced into the CNM phase, at least in the solar neighborhood. The numerical results (Figs. 
12a through 12c) show that this is true over most of the disk of the Galaxy. 
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6.4-2. Turbulently Supported WNM 



The discussion in the previous section is based on a highly idealized model of the ISM, 
in which the gas is supported entirely by thermal pressure. In fact, turbulence makes a 
substantial contribution to the support of the gas; cosmic rays and magnetic fields appear 
to make less of a contribution to support of H I near the plane since their scale heights are 
much greater than that of the CNM and the gradient therefore weak, though their substantial 
pressures must somehow be anchored by the weight of the ISM (Boulares & Cox 1990). The 
increase in the scale height reduces the mean density of the gas, and hence the inferred 
thermal pressure. The turbulence leads to large pressure fluctuations, and these will drive 
some of the gas into the cold phase (e.g., Hennebelle & Perault 1999, 2000). The condition 
for a substantial amount of CNM in a turbulent ISM is therefore that the thermal pressure 
exceed P m i n , which ensures that the gas that is compressed into the CNM can remain there. 

We analyze first the solar circle and compare the thermal pressure to P m i n assuming 
all H I is WNM. We redistribute the 2.75 x 10 20 H I atoms cm' 2 of CNM into the WNM 
components found by Dickey & Lockman (1990); the result is a gas with a central density 
of 0.31 cm~ 3 . This procedure should provide a lower bound on the thermal pressure of 
this hypothetical ISM in which the H I is pure WNM, since we are assuming that there is 
enough additional turbulent energy injection to lift the mass in CNM up to the height of the 
WNM. With this assumption, the resulting thermal pressure is P/k = 2700 K cm~ 3 . Since 
this significantly exceeds P m m/k = 1960 K cm -3 and since the local ISM is turbulent, we 
conclude that it must be in two phases. 

We can make analogous arguments to the rest of the Galaxy to show that, assuming all 
of the H I were WNM, the thermal pressure would exceed P m i n , violating our assumption, 
and that therefore CNM must exist. We begin by determining a lower limit on the thermal 
pressure produced by WNM gas at other positions in the Galaxy. In § 3.1, we have presented 
the mean density (nni) of H I in the plane rather directly derived from the observations. 
Assuming that the volume filling factor of WNM is much greater than that of the CNM, 
then 



where Mcnm/^wnm is evaluated in the midplane. The thermal pressure in the medium is 
given by 
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-Pth.WNM — l-l^WNM&^WNM • 

In terms of the "observable" (nni), -Pth.wNM can be written as 



(35) 
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D _ ^ ^ ("Hl)fcr wNM _ (-PwNm) ^fij 

/hi(1 + ^cnm/^wnm) /hi(1 + ^cnm/Mvnm) 



The "mean" thermal pressure (Pwnm), the thermal pressure assuming the denominator 
is unity (for example, the observed mean density is all WNM and /hi — 1), is shown in 
Figure 12d (assuming Twnm = 8000 K) along with P m i n , Pmax, and Pth.ave for our standard 
N c \ = 1.0 x 10 19 cm" 2 case. Following the distribution of H I in the Galaxy, (Pwnm) is flat 
with galactocentric radius out to 13 kpc where it starts to drop exponentially. 

As we did in the non-turbulent case, we can make the ansatz that all the H I gas 
is WNM and ask whether the thermal pressure exceeds P m i n or P max - In contrast to the 
argument above for the local ISM, we assume that there is no additional turbulent energy 
injection to raise the height of the converted CNM. In this case (Pwnm) is a lower limit to the 
thermal pressure since /hi may be less than unity (although it must be borne in mind that a 
hypothetical ISM in which the H I is pure WNM is likely to have / H i closer to unity than the 
actual ISM). Comparing (Pwnm) to the model thermal pressure curves, we conclude that in 
this case the thermal pressure would exceed P max for Galactic distances 8 kpc ^ R 16 kpc; 
our ansatz is violated; and a two-phase medium must exist in the outer Galaxy. Moreover, 
over much of the Galaxy (Pwnm) exceeds P min and hence pressure fluctuations in a turbulent 
ISM will produce a two-phase medium. 

We also note in Figure 12d that (Pwnm) falls everywhere below the non-turbulent pres- 
sure Pwnm' derived from the weight of the H I in the galactic potential. Essentially, this high 
value of Pwnm' derives from the assumption of no turbulence and therefore has a thermal 
scale hight of ~ 80 pc (see Appendix A), whereas (Pwnm) takes into account the observed 
half height (~ 265 pc) of the H I and reflects the importance of turbulence for the dynamics 
of the interstellar gas. 



7. An Analytic Model for Two-Phase Equilibria 

In order to understand the physical reasons for the results just presented, and to obtain 
scaling laws, it is convenient to have an approximate analytic model. The details of this 
model are presented in Appendix C; we summarize here our procedure and the results. 

The thermal balance pressure P is found as a function of T by equating a simple 
analytic equation for the total ([CII] +[OI]) cooling to the analytic expression for grain/PAH 
photoelectric heating. The simple cooling equation holds for 100 K < T < 1000 K, which 
marks the temperature range of validity of our analytic solutions. Several other parameters 
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enter the pressure equation, including the FUV field G' Q , the total gas ionization rate by 
cosmic rays and EUV/X rays Q, the dust/PAH abundance Z' d , and the gas phase metal 
abundance (especially of the coolants C and O) Z' g . The primes denote normalization to 
the local values given in Table 1, so that all these parameters have value unity at the solar 
circle. In other words, G' Q = Gq/1.7, and Q = Ct/10~ 16 s _1 , where Q includes primary and 
secondary ionizations. The grain photoelectric heating depends on the charge state of the 
PAHs/grains, which in turn depends on the electron density n e in the gas. At a typical 
column N c \ = 10 19 cm -2 , H + (not C + ) supplies the electrons. Therefore, n e ~ n(H + ) and we 
solve for n(H + ). The source of H + is the photoionization of H by EUV/X-rays, and the sink 
is recombination with negatively charged PAHs, or PAH~. PAH~ is produced by electron 
attachment on neutral PAHs (neutral PAHs dominate the PAH population in the parameter 
space valid for our analytic P expression), and is destroyed by photodetachment in the FUV 
field (G' Q ). 

With an analytical expression for P, we then take dPj dT = to find the temperature 
T( min ) at the pressure minimum, which substituted into our expression for P gives P m in- We 
put the parenthesis in the subscript for T( m i n ) to emphasize that this temperature is not a 
minimum temperature; in fact it is the maximum temperature of the CNM! The hydrogen 
density at P min is then n min = P min /[l.lA;T( min )]. Appendix C describes all the essential 
and competing reactions, and gives analytic expressions for the cooling, heating, PAH~ 
abundance, n e , and P. We summarize here the results for P m i n , T( min ), and n min , which 
we emphasize are valid over a range of parameter space centered on solar circle values, as 
discussed below. We find 

Q' ( -&\ 

P miQ /k = l.ln min T (min) ~ 8500 , ^oaes cm " 3 K > ( 37 ) 

l + 3.l(^J 

T (min) ~ 243 K , (38) 

and 

^mi„ * 31 °° j x 365 cm" 3 . (39) 

1 + 3.1 ' " z 



CI 

These equations are valid over the following range of parameters. The first condition is that 

4. 6 x 10-2 <^p< ii. (40) 

The lower limit is set so that the grain/PAH photoelectric heating is significantly affected 
by the positive charging of the grain/PAHs, which simplifies the analytic expression for the 
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gas heating. It also ensures that most PAHs be neutral and not PAH - , and that FUV 
photodetachment destroys PAH~ and not reaction with H + ; these assumptions were made 
in deriving the analytic expression. The upper limit is set so that H + is destroyed by PAH - , 
and not neutral PAHs. It also ensures the less restrictive condition that neutral PAHs, and 
not PAH + dominate the PAH population. The second condition is that 

/ r i\ 0.27 

6.5 x 1(T 3 4 2 < Z™<I? PAH <4.1. (41) 

The upper limit ensures that the gas-phase abundance of H + exceeds that of C + (H + supplies 
the electrons in the gas). As an interesting sidenote, we find that there is no P m i n in the 
temperature range where [C II] 158 /iin and [O I] 63/iin dominate the cooling, if the electrons 
are supplied by and if the grain heating is significantly affected by positive charge. The 
ionization of H and He is therefore generally critical to the two-phase phenomenon. The lower 
limit ensures that atomic hydrogen collisions, and not electrons, dominate the excitation of 
the [C II] line. The final condition is that 

\ 0.365 

£38*W (42) 

This condition assures that H + is destroyed by PAH~, and not by recombining with gas 
phase electrons. We found it interesting that this is generally the case, for the wide range of 
conditions centered on local values. 

Given that all the above conditions are met, equations (37) through (38) give not only 
the absolute values of P m i n and T( min ), but the scaling of P m i n with the parameters G' , Z' d , 
and Z' g . It might be noted that Z' d and Z' g may often scale linearly with respect to each other, 
and with the total elemental metallicity Z. It is also noteworthy that in this simple analytic 
solution, T( min ) is independent of all the parameters Z' , Z' d , G' , and It depends primarily 
on the gas temperature dependence of the cooling rates and the grain photoelectric heating 
rate. In much of parameter space, the numerical solution gives T( min ) ~ 180 — 350 K, which 
lie within about ±40% of the constant 243 K analytic result. The numerical calculation 
of T( min ) is sensitive to small contributions from processes we have ignored or simplified in 
the analytic treatment (e.g., electrons supplied by C + , our approximate cooling function, a 
simplification of the grain photoelectric heating). However, our analytic results for P min and 
the scaling of P min are much more robust; P min is not sensitive to the exact value of T( min ), 
since it lies at a minimum with respect to changes in T. 

Figure (13) illustrates the region of validity for the analytic solution for P min (eq. [37]). 
We plot the conditions given in equations (40) through (42) as functions of Z' = Z' g = Z' d and 
G' /Q. The constraint that hydrogen collisions dominate the excitation of [C II] (lower limit 
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in eq. [41]) is readily met and not a factor as long as the other conditions are satisfied. We 
have performed several checks of the analytic solution for P min and T( min ) against the detailed 
results of our numerical code. We held three of the parameters fixed at the local values (= 1 
in the notation of the analytic equations), and varied the fourth over a factor of 100 from 
0.1 to 10. This test showed that over this range the absolute value of the analytic P m i n was 
within ±50% of the numerical value as long as Z' g was less than 3 and Z' d was between 0.3 
and 5.0. The scaling of P min was accurate to about ±45%, and T( min ) was accurate to about 
a factor of 2.5. In the numerical runs, T( min ) varied from 180 to 630 K, with a value of 258 
K at the solar circle. This test, however, occasionally violated the conditions of validity. If 
we restrict the test strictly to the regime of validity, the agreement for P min is unchanged, 
but the range of T( min ) in our numerical calculation is reduced to 180-465 K, with much of 
parameter space 180 — 350 K. 

As another test of our analytic solution we compared P m i n to our numerical results in 
which we varied Z' g = Z' d from 0.01-10 and G' /Q from 0.1 to 300. The shaded region in 
Figure (13) shows the range in which the analytic solution agrees with the numerical results 
to within ±50%. It is clear that the applicable range of validity extends well beyond that 
given by the restrictive conditions expressed in equations (40) through (42). As a final test 
we compared the analytic solution for P min with the calculated P min as a function of Galactic 
radius (§ 6.1 and Table 3). We find that the calculated P min can be fit with an exponential 
in Galactic radius as 

P min = 1.1 x 10 4 exp(-P fc /4.9) K cm" 3 . (43) 

This fit is good to within ±17% between R k = 3 and Rk = 18. Substituting values for G' , 
(' t , Z' d , and Z' g from Table 2 into equation (37) yields analytical results which are good to 
within ±15% of the numerical results. 



8. Comparison with Observations and Previous Work 

8.1. Comparison with Local Observations 

In this section we compare our model results with three critical observations of the CNM: 
the thermal pressure, the gas temperature, and the C + cooling rate per hydrogen atom. We 
also briefly examine the constraints imposed by the observed C I/C II ratio. Our model 
results for the local Galaxy and the effects of various model parameters are summarized in 
Tables 3 and 4. The results presented in this paper so far have been for the WNM/CNM 
interface at a cloud column of N ci . We also list in Table 4 results for a CNM cloud interior 
at a depth of 1 x 10 20 cm -2 . The mean CNM column density measured by Heiles & Troland 
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(2002) is 5 x 10 20 cm -2 . For a slab cloud, a typical point is at a depth of 1/4 the observed 
column density (~ 1.25 x 10 20 cm -2 ), where half the mass is at greater column and half at 
lower column. For a spherical cloud, the average line of sight passes through a column of 
2N/3, where N is the column density through the cloud center, and the half mass depth is 
at 0.1 iV from the surface. This half mass depth corresponds to ~ 0.75 x 10 20 cm -2 . Our 
fiducial column is between the two limits for slab and spherical clouds. 

The thermal pressure in the ISM has been measured through studies of the population 
distribution of the C I fine-structure levels based upon FUV absorption lines originating 
from these levels (Jenkins et al. 1983; Jenkins & Tripp 2001). By necessity, such studies are 
largely limited to the local solar neighborhood. Jenkins & Tripp (2001) find a mean observed 
C I*/C I to t ratio of 0.196 and derive a mean pressure of P th = 2240 K cm -3 from this ratio. 
This pressure is somewhat lower than the previous result of Jenkins et al. (1983) who found 
Pth ~ 4000 K cm" 3 . Our standard parameter set yields P min = 1960, K cm" 3 P max = 4810 K 
cm" 3 and Pth, avc = 3070 K cm" 3 and our average pressure would appear to be higher than the 
observed value. However, Jenkins & Tripp (2001) note that the derived pressure is sensitive 
to the CNM temperature and atomic constants where they used a CNM temperature of T 
= 40 K and radiative decay rates from Froese Fischer & Saha (1985). The more recent 
radiative rates of Galavis, Mendoza, & Zeippen (1997) are approximately 10% higher than 
those used by Jenkins & Tripp (2001) and our model temperature of 71 K in the CNM cloud 
interior (see Table 4) is considerably higher than 40 K. Both of the these differences tend to 
increase the derived thermal pressure for the same observed population ratio. (We note that 
the excitation rates of C I by collisions with H I from Launay & Roueff 1977 are expected to 
be accurate to within 10%; E. Roueff 2002 private communication). Including the updated 
radiative rates, plus the optical/UV pumping of the C I lines as discussed in Jenkins & Tripp 
(2001) and Jenkins & Shaya (1979) we find a C I*/C I to t population ratio of 0.190 at the 
CNM cloud surface, and a C I*/C I to t ratio of 0.201 in the cloud interior (Table 4). These 
values are within 3 percent of the observed ratio of 0.196. The "Low 0pah" (0pah = 0.25) 
model compares least favorably to the observations with a population ratio of between 0.145 
and 0.162. Although the differences between models are not large, we conclude that our 
standard model with an average thermal pressure of P t h,ave = 3070 K cm" 3 provides the 
best agreement with the observed C I*/C I to t ratio. 

In addition to the thermal pressure, observational tests are provided by the CNM tem- 
perature and C + cooling rate. Heiles (2001), using H I emission/absorption experiments 
along 19 lines of sight, found that of the total CNM column detected, most (~ 61%) lies 
in a narrow temperature range between 25 K and 75 K with a peak near 50 K. This value 
was recently revised by Heiles & Troland (2002) when an additional 60 lines of sight became 
available and it was possible to derive temperatures for separate low latitude < 10°) and 
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high latitude > 10°) directions. They find a mean mass weighted (low latitude) CNM 
temperature of 99 K and a median temperature of 63 K with half the column density above 
and half below this temperature. 

Additional measures of the gas temperature of diffuse clouds come from UV absorption 
line observations of the J = and J = 1 level populations of H 2 in the ground vibrational 
state. The excitation temperature of these levels (sometimes denoted T 01 ) will equal the 
gas kinetic temperature as long as collisions with H + are sufficiently rapid to thermalize the 
H 2 ortho/para ratio, and that other processes (e.g., FUV pumping or H 2 formation) do not 
drive the ratio away from thermalization (see e.g., Dalgarno, Black, & Weisheit 1973; Black 
& van Dishoeck 1987; Burton, Hollenbach, & Tielens 1992; Sternberg & Neufeld 1999). We 
estimate that the ionization provided by the X-ray/EUV field provides sufficient H + so that 
T 01 is a good measure of the gas kinetic temperature in CNM clouds. Early estimates of the 
gas temperature from Copernicus observations found T 01 = 77 ± 17 K (Savage et al. 1977) 
while recent estimates using FUSE report Tqi = 68 ± 15 K (Rachford, et al. 2002). Of the 
models shown in Table 4 the "High 0pah" (0pah = 1) calculation has temperatures in the 
cloud interior more than 1 sigma higher than the H 2 temperatures. Our CNM temperature 
of 71 K for the standard model lies between the median and mean H I temperatures of Heiles 
& Troland (2002) and is in good agreement with the H 2 measurements. 

We note that Heiles & Troland (2002) found ~ 4% of the CNM column to lie at tem- 
peratures less than 25 K. Such extremely cold H I has been also reported, for example, by 
Gibson et al. (2000) who find local H I gas at T = 16 - 32 K, and by Knee & Brunt (2001) 
who find an H I supershell at it! = 16 kpc with T = 10 K. Our standard model is unable 
to produce these low CNM temperatures. Such low temperatures might be produced in 
gas depleted in PAHs, which reduces the photoelectric heating (as suggested by Heiles & 
Troland 2002), or in gas that contains a sufficient molecular H 2 abundance to cool to these 
temperatures. 

The reported [C II] 158 /im cooling rate per hydrogen atom varies widely depending on 
the method and direction of observation. The two main methods are derived from (1) UV 
absorption line measurements of the column density of C II*, N(C IP), where C IP means 
C+ 3 P 3/2 , and (2) IR observations of the [C II] 158 urn line intensity J([C II]). From JV(C IP) 
along a line of sight, and the total (H I + H II + 2H 2 ) column of hydrogen nuclei, iV(H), the 
cooling rate is found from nA = E21A21N '(C IF)/iV(H) where E21 and A21 are the energy and 
radiative decay rate for the [C II] 158 /zm transition. The total hydrogen column is usually 
inferred from the column of S + . For the IR method, the results are usually reported using 
H I 21 cm emission to obtain the column of neutral hydrogen iV(H I). From the intensity 
of the [C II] line and iV(H I), the average cooling rate is given by nA = 47r/([C II])/7V(H I). 
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Note that the UV and IR measures are not directly comparable since the former is per H 
nucleus and latter is per H I atom. One source of uncertainty in both methods is the amount 
of WNM gas along a given line of sight. This is because the total column of hydrogen arises 
from both WNM and CNM gas whereas the C II* resides in mainly the CNM gas alone. 
Thus the effect of WNM gas is to lead to an underestimate of the derived cooling rates in 
the CNM. This problem is especially severe for the high latitude IR observations in which the 
H I column extends several hundred parsecs above the plane where there is not much cold 
gas. There are also processes which may lead to overestimates of the derived cooling rate. 
The UV absorption method can be influenced by radiation from the background stars which 
provides a UV radiation field, and subsequent gas heating rate, greater than the average 
interstellar field. The IR line and UV absorption methods can be contaminated by WIM 
gas, which has a cooling rate per hydrogen in the [C II] transition similar to the CNM so that 
the WIM contributes to both IR emission and UV absorption. (However, generally along a 
line of sight the column of WIM is much less that that of CNM). 

The UV absorption method was used in the Galactic plane by Pottasch, Wesselius, & 
van Duinen (1979), who found an average value of nA w 1 x 10~ 25 erg s" 1 H -1 , and by 
Gry, Lequeux, & Boulanger (1992) who reported nA = 3.5^21 x 10~ 26 erg s _1 EU 1 . The 
Pottasch et al. (1979) value is probably biased towards higher rates due to UV illumination 
from background stars. On the other hand, the Gry et al. (1992) results are probably 
biased towards lower rates. Their lines of sight were selected to have low column density 
and are thus expected to contain more WNM gas on average then a typical line of sight. 
The C II* columns have been determined for a few high latitude lines of sight from which 
the cooling rate can be derived. For example, from Fitzpatrick & Spitzer (1997), we find 
nA = 0.8 — 3.5 x 10~ 26 erg s _1 H _1 in the cool components toward HD 215733, and Savage 
et al. (1993) find nA = 1.4 x 10~ 26 erg s" 1 FT 1 toward 3C 273. The Savage et al. (1993) 
observation is probably a lower limit to the CNM cooling rate. The direction towards 3C 
273 is a well studied line of sight with low mean density and much of the C II* may reside in 
WNM gas. In addition, the hydrogen column (inferred from S + ) samples the entire column 
throughout the disk and halo and thus overestimates the column in CNM alone. 

The IR line emission method has been applied to high latitude directions. Bock et al. 
(1993) and Matsuhara et al. (1997) using a balloon borne experiment measured the 158 jum 
line emission in a path covering Galactic latitudes b = 33° to 50°. Matsuhara et al. (1997) 
obtain a cooling rate of 1.6 ± 0.4 x 10~ 26 erg s _1 (H I) -1 (i.e., per neutral hydrogen) for 
columns A^(H I) < 2 x 10 20 cm~ 2 . Bock et al (1993) using the same C II observations, 
but for V(H I) > 10 20 cm" 2 and excluding anomalously low C II/continuum ratios, found 
2.6 ±0.6 x 10~ 26 erg s _1 (H I) 1 . These rates are slightly lower than the COBE observations 
of 2.65 x 10~ 26 erg s" 1 (H I) -1 (Bennett et al. 1994). 
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We can estimate the correction to these high latitude IR observations due to WNM and 
WIM gas (see discussion in WHMTB). The C + line intensity arising in the ionized gas can 
be estimated from the Ha line intensity available from the WHAM data set (Reynolds et 
al. 1999) 3 . We find that the Ha line intensity varies from about 0.5 to 1 R over the path 
of the Bock et al (1993) experiment. Using /(Ha) = 1 R and the conversion between /(Ha) 
and /([C II]) (e.g., Reynolds 1992) we find /([C II]) = 1.5 x 10" 7 erg s" 1 cm" 2 sr" 1 from the 
ionized gas. Subtracting this from the Matsuhara et al. (1997) cooling rate we get roughly 
1.1 x 10" 26 erg s" 1 (H I)" 1 . Following WHMTB we use the average WNM and CNM columns 
derived from the observations of Kulkarni & Fich (1985) yielding iV" C NM = 0.4 x 10 20 esc \b\ 
cm" 2 and N WNM = 1.0 x 10 20 esc |6| cm" 2 . From WHMTB equation (12) and neglecting 
the neutral component of the WIM we find a cooling rate in CNM of 3.9 x 10" 26 erg s" 1 
H" 1 . We conclude that although the rates are quite uncertain due to contributions from 
WNM and WIM gas, a conservative lower limit is that the [C II] cooling rate in the CNM is 
greater than ^ 3.0 x 10" 26 erg s" 1 H" 1 . All of the models presented in Table 4 satisfy this 
constraint although the low 0pah, low PAH abundance, and low Go cases barely do so. 

We finally consider the model constraints imposed by the C I/C II ratio. From the 
C I and H I data compiled in Welty & Hobbs (2001), and assuming a constant C II/H I 
ratio of 1.4 x 10~ 4 , we find that C I/C II ~ 3 x 10~ 3 . The abundance of neutral metals, 
as typified by the C I/C II ratio, strongly rises with increasing 0pah through the direct 
recombination of C + on PAH", and the C I/C II ratio increases by a factor 8.2 as 0pah 
increases from 0.25 to 1. The effect of a lower PAH abundance is to diminish both the ion 
recombination on grains and the photoelectric heating rate, resulting in both lower thermal 
pressures and lower C I/C II ratios. Decreasing the FUV field from 1.7 to 1.1 decreases 
the rate of photoionization of C I and the C I/C II ratio rises by a factor of 1.6. We find 
that the "High 0pah and "Low Go" models can be safely ruled out by the required C I/C II 
ratio. We shall present in a future paper a more detailed analysis of the implications of the 
C I/C II ratio for translucent cloud models. In summary, we find that the average thermal 
pressures are not dramatically affected by the range of parameters shown in Table 4, however 
in a detailed comparison with observations we find (1) "Low 0pah" models poorly match 
the C r/C tot ratio and C + cooling rate, with values lower than observed (2) "High 0pah" 
models produce temperatures and C I/C II ratios higher than observed, (3) "Low PAH" 
models produce low values for the C I*/C to t ratios and C + cooling rate, and (4) "Low Go" 
models produce low vaules for the C I*/C to t ratios and C + cooling rate, and high values for 
the C I/C II ratio. 



3 http://www.astro.wisc.edu/wham/index.html 
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The discussion in this section highlights the need for more observational and labora- 
tory study in the areas of PAH chemistry and thermal processes in the diffuse interstellar 
medium. In the near future several ground-based, airborne, and space based telescopes 
(ALMA, SOFIA, Herschel) may enable us to measure C I fine-structure level populations 
throughout the Milky Way as well as in other galaxies. These observations will provide 
stringent tests for our model and for the general importance of the thermal instability for 
the two-phase medium of interstellar gas. We note that SOFIA and Herschel will also allow 
us to study the [C II] 158 /im emission along specific sight lines in the Milky Way at high 
spatial and spectral resolution. Together, this will provide a powerful way to study the 
physical conditions in the ISM and the processes that dominate its structure. 



8.2. Comparison with Extragalactic Observations 

In comparison to the Braun (1997) observations, we find that based on a phase diagram 
analysis alone we do not expect a strict cut-off in the cold gas at the Galactic i? 2 5 radius 
(i?25 = 12.25 kpc). We also find a CNM temperature that decreases in the outer Galaxy, 
where as Braun (1997) finds peak brightness temperatures to increase. To reconcile these 
differences, we suggest that Braun (1997) is mainly detecting the photodissociated outer 
envelopes of molecular clouds. The inferred high column densities (N ^ 3 x 10 21 cm' 2 ) and 
temperatures (T ~ 80—150 K) reported by Braun (1997) are more typical of photodissociated 
gas near star forming regions than of the diffuse CNM clouds illuminated by the general 
interstellar radiation field on which we have been concentrating. A similar interpretation of 
the H I emission in M101 has been proposed by Smith et al. (2000). 

Using the photodissociation region models of Kaufman et al. (1999), modified by the 
chemistry and thermal processes discussed in this paper, we find that if the gas density and 
UV field illumination are fixed, the H I brightness temperature increases if the metallicity 
(and therefore the gas phase coolants and grain abundances) decreases as expected in the 
outer portions of galaxies. Near regions of star formation, the UV radiation field does not 
drop with the disk radial scale length, but remains relatively constant and depends mainly 
on the illumination from the nearest star-forming region. For local ISM metallicities, the H I 
21 cm line is only marginally opaque, so that the maximum brightness temperature is set by 
the cool interior regions of H I clouds. However, at lower metallicities (and lower dust-to-gas 
ratios) the same visual optical depth is achieved at larger H I columns and thus the H I 
to H2 transition is pushed deeper into the cloud. The larger H I column forces the 21 cm 
transition to become optically thick while still within the warm regions, thereby increasing 
the H I brightness temperature. For metallicities as low as 0.2-0.3 times cosmic, typical of 
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those at the P 2 5 radius in M101 and other spirals (Garnett, et al. 1997) we find that the 
temperature increases by 40 — 60 K, comparable to that seen by Braun (1997). 

As discussed by Sellwood & Balbus (1999), a fairly uniform H I 21 cm line width of ~ 6 
km s" 1 is seen in the outer disk of spiral galaxies. Instead of indicating WNM gas, Sellwood 
& Balbus (1999) suggest that the dispersion arises from a cold, MHD-driven turbulent gas. 
Based on the results of WHMTB, they argue that a 5000 K gas inferred from the line width, 
would be thermally unstable. Our temperature results for WNM in the the outer Galaxy 
shown in Table 3 confirm that 5000 K gas is indeed thermally unstable there. Furthermore, 
our thermal pressure arguments for cold gas in the outer Galaxy, are also consistent with 
their turbulence hypothesis. However, we also anticipate that the H I in the outer Galaxy is 
a two-phase medium with a substantial amount of WNM. 



8.3. Comparison with Previous Work 

Our results differ from those reported in Elmegreen & Parravano (1994), who found that 
CNM gas extends to only ~ 12 — 14 kpc in the Galactic plane. Using the thermal processes 
and chemistry from Parravano (1987), they estimated the range of thermal pressures (P m m < 
P < Pnax) over which two phases are possible, as a function of FUV radiation field and 
metallicity. They also estimated the midplane thermal pressure and concluded that the 
thermal pressure would drop below P min by a Galactic radius of ~ 12 — 14 kpc. 

The difference between our results and theirs arises in large part due to the approxima- 
tions they made in estimating the pressure. Following Elmegreen (1989), they wrote the total 
pressure, P tot , as being proportional to the product of the gas surface density £ gas and the 
total (gas plus stars) surface density, £ gas+sta rs, or P tot oc £ gas E gas+stars . They then assumed 
that the total surface density is proportional to the gas surface density, S gas+Stars oc £ gas , so 
that P tot oc S gas , and that the ratio of thermal to total pressure is constant so that P t h oc S gas . 
Finally, and most important, they assumed that S gas oc exp(— P/PTr) with Hr ~ 4 kpc. As 
a result, they obtained Pth/k ~ 10 5 ' exp(— 2R/Hr), which results in a thermal pressure of 
only P/k ~ 30 K cm -3 at 16 kpc. In fact, as discussed in § 3, the H I surface density remains 
above 5 M pc~ 2 out to R ~ 15 kpc. Our value for the pressure of a pure WNM, Pwnm', 
which is based on the observed H I surface density, is a factor 67 higher than theirs at 16 
kpc. This higher pressure extends the region of cold gas to much greater radial distances 
than found by Elmegreen & Parravano (1994). 
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9. Summary and Discussion 

9.1. Model Assumptions 

Our model for the gas heating and ionization in the Galactic disk is based on a cosmic- 
ray rate that is constrained by comparisons of observations to chemical models (e.g., van 
Dishoeck & Black 1986; Federman, Weber, & Lambert 1996; van der Tak & van Dishoeck 
2000), direct observations of the local FUV field and soft X-ray field, and a theoretical 
estimate of the EUV intensity. The soft X-ray and EUV intensity is partly derived from the 
calculations of Slavin et al. (2000), who find that radiation from cooling supernova remnants 
can produce the fractional ionization seen in clouds at high latitude. The intensity used in 
this paper is an extension of their calculation to the Galactic plane. In addition, a stellar 
EUV component from Slavin et al. (1998) is added so that the total ionizing photon flux 
from the Galactic disk matches the recombination rate derived from Ha observations. An 
important component of the EUV and soft X-ray radiation transfer is the opacity produced 
by the WNM gas component. We have assumed that the opacity is provided by an ensemble 
of WNM clouds each of column density 7V cl (H I). The radiation fields at the WNM/CNM 
interface and cloud interior are found by passing the incident radiation through additional 
columns of iV c i(H I), and 1 x 10 20 cm -2 , respectively. 

Our confidence in the adopted local parameters is strengthened by our successful mod- 
eling of the densities and temperatures in the local WNM and CNM gas. Furthermore, we 
obtain a good fit to the [C II] cooling rate per hydrogen as derived by UV line absorption 
and IR line emission studies. The good match between theory and observations further in- 
dicates that we have included all the relevant heating and cooling terms in our model, an 
argument that is especially strong for the CNM phase in which [C II] dominates the cooling; 
any additional heating sources in our model would emerge as excess [C II] emission. This 
case is not as strong for the WNM in which the [C II] emission contributes only ~ 14% of the 
cooling and is a factor 16 weaker per hydrogen than in CNM. Additional heating terms do 
not strongly affect the fine-structure line emission due to the strong temperature regulation 
by Lja cooling. 

To obtain the FUV field in other regions of the Galaxy we have relied on observations 
of the gas surface density, metallicity, and OB star distribution, along with a numerical 
integration for the mean intensity. To check our assumptions, in § 6.3, we compared the 
infrared emission produced by dust grains heated by the interstellar radiation field with 
observations taken by the COBE satellite. The processes that determine the X-ray and 
cosmic ray distributions are certainly more complicated than our models allow. For example, 
the soft X-ray flux depends on the temperature and emission measure in the hot ionized gas 
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component and the optical depth to the emission regions. The temperatures and emission 
measures depend on the metallicity, and the optical depth depends on the structure of the 
interstellar medium. We note, however, that Slavin et al. (2000) find that the intensity of 
X-rays produced by SNR emission does not depend sensitively on the ambient density of the 
preshock gas. We have carried through our analysis by adopting a simple approach in which 
we use plausible arguments for scaling the soft X-ray and cosmic ray rates to other regions 
of the Galaxy based on the distribution of production sources (OB stars) and destruction 
sinks (various gas components). 

Note that we have chosen to extend the OB star distribution with a constant scale length 
Hft B = 3.5 kpc out to R — 18 kpc. Had we adopted the OB star distribution of Bronfman 
et al. (2000), P min and P max would have been lower and it would have been easier to form 
CNM in the outer Galaxy. Since the actual OB star distribution in the outer Galaxy is if 
anything below the one we have adopted, our conclusion that CNM must exist in the outer 
Galaxy is strengthened. 

Although the distribution of total H I surface density is constrained by observations, 
the separate column densities of CNM and WNM gas are not well determined away from 
the solar neighborhood. In calculating the opacity for EUV and soft X-ray photons, we have 
assumed that the ratio of WNM to CNM surface densities and scale heights are the solar 
neighborhood values (Dickey & Lockman 1990) throughout the Galaxy. This in turn implies 
that the ratio of CNM to WNM volume filling factors in the midplane is held constant with 
Galactic radius. A somewhat different prescription is given by Ferriere (1998) in which the 
volume fraction of WNM increases in the outer Galaxy. Note that if we allowed the WNM 
fraction to increase in the outer Galaxy, then the opacity to soft X-ray and EUV radiation 
would increase as well, thereby reducing P min and P max and making conditions less favorable 
for the existence of WNM. Another difference in our H I distributions is that Ferriere (1998) 
used a constant surface density of S H i = 5 M pc~ 2 in the outer Galaxy to R = 20 kpc where 
our distribution drops below 5 M pc~ 2 beyond R > 15 kpc. For constant surface density, 
at R = 18 kpc the pressure Pwnm' is a factor ~ 2 higher than for our H I distribution and 
would make it more likely that CNM gas can exist. In a future paper we shall compare the 
calculated [C II] emission with the observational data in an effort to independently derive 
the volume fractions of CNM and WNM gas. This will be particularly telling for the outer 
Galaxy, where the relative absence of photodissociation regions and H II regions should 
permit a clean distinction between WNM and CNM without the confusion of predominantly 
molecular or ionized gas. 
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9.2. Galactic Distribution of Two-Phase ISM 

As discussed in § 6.4, in Figures 12a through 12c we plot P m i n and P max as a function of 
position in the Galactic midplane for 3 values of the WNM cloud column N c \. We also plot 
the thermal pressure Pwnm' in the Galactic midplane that would result if all of the H I layer 
were WNM gas supported by thermal pressure. For the case in which there is no turbulence 
and the thermal pressure dominates the pressure, regions in which Pwnm' > -Pmax must have 
CNM gas. This is because only CNM gas can exist at these pressures, and mass will be 
converted from the WNM phase to the CNM until the pressure drops below P max , where a 
two-phase medium can exist. Our figures show that this condition is satisfied over much of 
the Galactic disk. 

Pwnm' is calculated assuming all the diffuse gas is WNM and that thermal pressure 
dominates and determines the vertical scale height, which we calculate in this case locally 
to be ~ 80 pc. However, the observed half-height of the 'WNM component" seen by Dickey 
& Lockman (1990) is ~ 265 pc. This result demonstrates that nonthermal pressure (due 
to turbulent motions, magnetic fields, and cosmic rays) dominates. We take an analogous 
approach to analyze the turbulent case. Assuming again that all the H I gas is WNM, 
we compare the thermal pressure to P m i n and P max - We use the observed (nm) (which 
includes the effects of turbulence and the greater scale heights which lowers (rim)) to estimate 
a lower limit (Pwnm) on the thermal pressure P t h,wNM = (Pwnm) //hi where /hi is the 
volume filling factor of the WNM (the rest is HIM). Since (Pwnm) exceeds P max in the outer 
(8 kpc ^ R ^ 16 kpc) Galaxy, CNM must exist in these regions. Since (Pwnm) exceeds P m i n 
(and turbulence likely drives the local pressures above P max occasionally), and since /hi < 1, 
we conclude CNM very likely exists at 3 kpc ^ R ^ 18 kpc. 

It is difficult, however, from our theoretical models, to rule out an interstellar medium 
with only HIM and CNM (and no WNM) in which the intercloud medium is filled with HIM 
that maintains a pressure P > P max on the CNM clouds. However, in such a scenario, the 
volume filling factor of the HIM must be nearly unity. If the HIM does not fill the intercloud 
medium, CNM would partially convert to WNM to fill the vacuum, and the pressure in the 
pervasive WNM would adjust such that P max > P > P m i n (see Parravano et al. 2002 for 
further discussion). We also note that the observation of H I 21 cm emission and absorption 
throughout the Galaxy strongly suggests the presence of a pervasive WNM. 

Assuming that P lies between P min and P max , we can use our models to predict the 
average thermal pressure from equation (26), and the results are given in Table 3. An 
approximate analytic fit to these results for our fiducial column density of N c \ = 10 19 cm~ 2 
is 

Pth,avcA = 1.4 x 10 4 exp(-P fe /5.5) K cm- 3 . (44) 
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We find that at fixed Galactic radius, the pressure does not change by more that a factor 
~ 3 over our range of cloud columns. For our fiducial column, the pressure drops from about 
8,200 K cm -3 at 3 kpc to 3100 K cm" 3 at the solar circle and to 600 K cm" 3 at 18 kpc. The 
drop in the thermal pressure from 3 kpc to the solar circle (a factor 2.7) closely matches the 
drop in the magnetic pressure inferred from radio observations: Beck (2001) estimates that 
B drops from about 10 /xG at 3 kpc to 6 fiG locally, corresponding to a pressure drop by a 
factor 2.8. (Note that these values for the field are larger than the rms field; as Beck points 
out, his values of B are (P 3 - 9 ) 1/3 - 9 .) 

We can now test the validity of our assumption that the turbulence parameter T = 
^cooiAshk ~ 1, so that it is meaningful to discuss a two-phase medium. Our results show that 
locally T « 0.1 for CNM and T « 0.3 for WNM at a pressure of P th /k = 3000 K cm" 3 
(see eq. [5]). As a function of Galactic radius we find that t coo i oc T/(nA) oc exp(P fc /2.94). 
Ignoring the weak variation of t s hk due to the variation in Swnm and in n^^, we have 
^shk oc <r<T^ oc exp(P fc /3.5). As a result, we find T oc exp(P fc /18.4), and we conclude that the 
turbulence parameter is ^ 1 and weakly dependent on radius throughout the Galactic disk. 

In Appendix B we discuss the potential role of turbulence in heating the WNM and CNM 
phases. Using the admittedly uncertain turbulent heating rate as a function of Galactic 
radius given by equation (B5), at R — 17 kpc we find that P m i n , Pth,ave, and P max are 
about a factor 2 higher than for the non-turbulent heating case. However, since (Pwnm) 
remains above P m i n to R — 18 kpc, turbulent heating does not change our conclusion that 
turbulent fluctuations will produce cold gas that is thermally stable in the outer Galaxy. The 
rate of turbulent heating does not exceed the rate of photoelectric heating out to R ~ 18 
kpc, and from equation (B3) we conclude that our assumption of thermal balance remains 
approximately valid. 

We conclude by summarizing our most important results. We have shown that both 
observational evidence and our theoretical models presented here indicate that the thermal 
pressure in the ISM of the Galaxy lies in the relatively narrow range between P min and P max 
for 3 kpc < R < 18 kpc. We have calculated P min (P), P max (P) and an estimate of the 
thermal pressure P t h,ave(P) in the Galaxy. We have shown that CNM gas must exist out 
to 18 kpc. We present phase diagrams for several galactocentric radii and for several cases 
of varying opacity to EUV and soft X-ray flux. Understanding the neutral phases of the 
ISM and their dependence on the radiation field is an important step in understanding the 
formation of molecular clouds and the global star formation rates in a galaxy. 
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A. Analytic Solution for Thermally Supported WNM at Solar Circle 

We can treat analytically the problem discussed in § 6.4.1 for the thermal support of 
the WNM if we assume that the mass in the disk is distributed exponentially with height 
above the plane, 

Pt = Ptoe- z/zd , (Al) 

where Zd is the scale height and p t is the mass density of all the matter — stars, gas, and dark 
matter. The dominant contribution to the mass in the disk at R = Rq is the "thin disk" 
of stars with z d = 180 pc (Dehnen & Binney 1998). Solution of the equation of hydrostatic 
equilibrium for the gas, which is assumed to be isothermal with sound speed a t h, gives 

n = n exp [-(3x + (3(1 - e~ x )] , (A2) 

where x = z/zj, and 

°th ™ 

In the limit of large x (z >> z d ~ 180 pc), the distribution becomes an exponential with 
a scale height z d j (3 = h 2 / 'z d ; in the limit of small x (which is more relevant for the gas 
distribution in the disk), the gas density approaches a Gaussian, 

n^n exp^-^^. (A4) 

Numerically, we have h = 0.89T 1 / 2 pc at the solar circle, based on a stellar density in the mid- 
plane of 0.115 M & pc -3 from Dehnen & Binney (1998) (with the 2.5 M & pc~ 2 augmentation 
to the thin stellar disk discussed above) and a gas density in the midplane of 0.034 M© pc~ 3 
from Dickey & Lockman (1990) and McKee (1990). For an adopted temperature of 8000 K, 
this gives h = 80 pc for the WNM. 

As remarked above, the value of the midplane density is determined by requiring that 
the column density agree with the observed value, 

POO 

N = 2n hp 1/2 / exp [-fix + 13(1- e^)] dx = (2n) 1/2 n h(f) /3 . (A5) 
Jo 

The fact that the mass is distributed exponentially rather than uniformly increases the 
Gaussian scale height of the gas from h to <f>ph. A little algebra shows that the factor (pp is 
approximately 

1 / 2 \" ! 1 



1+ sUJ + m (A6) 
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This approximation is accurate to within 2% for (3 > 1. For /3 > 1 we have 1 < 0^ < 1.35, 
so the deviation from Gaussian behavior is not large. With these results, we then find that 
the pressure the H I would exert in the midplane if it were all WNM is 



For our adopted parameters (JV H i = 6.2 x 10 20 cm" 3 , z d = 180 pc, T = 8000 K and h = 80 
pc), we find -Pwnm'A = 7800 K cm" 3 at the solar circle. By comparison, the exact numerical 
solution using the Dehnen & Binney (1998) potential gives -Pwnm'A = 8615 K cm" 3 , which 
is satisfactory agreement. 



In addition to the heating processes we have considered — photoelectric, X-ray, and 
cosmic-ray heating — turbulent, or mechanical, heating may also be important. Cox (1979) 
estimated that about 30% of the energy of a supernova would go into the compression of 
interstellar clouds followed by radiative losses. Numerical simulations (Cowie, McKee, & 
Ostriker 1981) confirmed this estimate. Spitzer (1982) estimated that about 4% of the en- 
ergy of a supernova would go into acoustic waves, and suggested that the absorption of 
these sound waves could be an important source of heating for the warm phase of the ISM. 
Ferriere et al. (1988) generalized this discussion to consider the generation and damping of 
hydromagnetic waves produced by supernova remnants (SNRs). Minter & Spangler (1997) 
used observations of interstellar scintillation to infer the amplitude of fluctuations in the 
ISM, and then set constraints on how these fluctuations damp. Minter & Balser (1997) and 
Mathis (2000) studied the effect of a turbulent heating rate of T ~ 10" 25 erg s" 1 H" 1 on 
the WIM. Sellwood & Balbus (1999) have estimated the heating rate due to the dissipation 
of turbulence generated by Galactic differential rotation. Turbulent heating has also been 
investigated in molecular clouds (e.g., Stone, Ostriker, & Gammie 1998; Mac Low 1999). 

Despite more than two decades of work on turbulent heating of the ISM, the rate 
remains very uncertain. This uncertainty stems directly from our lack of understanding of 
interstellar turbulence — how it is generated, how it propagates, and how it dissipates. We 
have already encountered this uncertainty when we tried to estimate T, which measures the 
degree to which non-turbulent heating is in balance with radiative cooling (see § 2). Most 
of the heating associated with SNRs is very intermittent, with gas being heated by shocks 
that are separated by long time intervals. As a result, even though about 30% of the energy 
of a SNR may go into cloud heating, most of this heat may be radiated away in a short 
time while the gas is substantially hotter than average. This is particularly true for the 
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B. Turbulent Heating 
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WNM, which has a rapidly rising cooling rate above 10 4 K due to Lya cooling. The fact 
that Heiles & Troland (2002) did not find many WNM features with line widths above that 
corresponding to the temperature we have calculated for the WNM suggests that impulsive 
heating is not a dominant process in determining the temperature of the interstellar H I. Cox 
(1979) reached the same conclusion based on the fact that the observed level of turbulence 
in clouds is relatively small. 

In view of the uncertainties associated with turbulent heating in the ISM, we have not 
included it in our basic models. Here we shall estimate the rate of turbulent heating and 
determine how it would affect our results. 



B.l. Dissipation of Turbulent Energy 

Let e be the rate of dissipation of turbulent energy per unit mass, which is equivalent 
to the turbulent heating rate per unit mass. On dimensional grounds, we expect e ~ Sv 3 /£, 
where 8v is the rms velocity in a region of size £ (Landau & Lifshitz 1987; Stone et al. 1998; 
Mac Low 1999). In other words, the kinetic energy per unit mass is dissipated in a time of 
order £/5v. A simple global estimate shows that turbulent heating is in fact unimportant 
in the ISM: The turbulent energy in the ISM is about O.bMSv 2 ~ 10 54 erg, where we set 
M = 10 9 M Q and 5v = 10 km s _1 . The scale on which the turbulence in the WNM has an 
amplitude of 10 km s" 1 is about 200 pc (§ 2), so the turbulent dissipation time is about 20 
Myr. 4 The resulting heating rate is about 1.6 x 10 39 erg s _1 , or 4 x 10 5 L Q . By comparison, 
the luminosity of the Galaxy in the C II 158 /im line is 5 x 10 7 L & (Wright et al. 1991), so 
turbulent heating is negligible on a galactic scale. However, as we shall see, it is relatively 
more important in the outer Galaxy. 

We make the dimensional argument for the heating rate exact by introducing the con- 
stant <f>£, 

= 3 3/ ^ j , (Bl) 

where a is the ID turbulent velocity dispersion. The turbulent heating rate per hydrogen, 
r turb , is then given by 

3 

Tturb = - = fine = 3 3 / °— , (B2) 



4 At first glance, it is puzzling that this estimate for the dissipation time exceeds the estimated time 
interval between shocks in the WNM of about 5 Myr in §2. However, as discussed above, much of the shock 
energy is radiated promptly and does not contribute to the level of subsonic turbulence. 
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where nu = 2.34 x 10 24 g is the mass per hydrogen. 

The numerical simulations of Stone et al. (1998) and of Mac Low (1999) show that 
equation (Bl) applies to both subsonic and supersonic turbulence. In supersonic turbulence, 
the energy is dissipated primarily in shocks, so the dissipation is highly localized in space 
and time. As discussed above, for this reason shocks do not contribute effectively to the 
general heating of the CNM and WNM, so we are more interested in the subsonic case. 
In the subsonic, non-magnetic case, we expect the turbulence to follow the Kolmogorov 
scaling, in which e is independent of scale (a oc i 1 / 3 ). When magnetic fields are included, 
Goldreich & Sridhar (1995) found the same scaling, so this should be generally valid in the 
ISM for scales such that the flow is subsonic. Since e is independent of scale for subsonic 
turbulence, the turbulent heating in this case should be widely distributed, albeit with 
substantial fluctuations associated with intermittency. On the other hand, for supersonic 
turbulence, we expect a oc £ q with q > 1/3 (Larson 1979; Boldyrev 2002). This scaling is 
observed in molecular clouds (Heyer & Schloerb 1997), which are highly supersonic. In this 
case we have e oc a oc Ai 3 ^ 1 ^, where Ai is the Mach number of the flow. The fraction of the 
total turbulent heating that is subsonic and thus widely distributed is thus about l/Ai 3 ^ 1 ^ 
for Ai > 1, corresponding to 1/A4 for q — | (Larson 1979) and 1/Ai 1 ^ 3 for q = 0.375 
(Boldyrev 2002). The length scale that separates subsonic from supersonic turbulent motions 
in the WNM, and of CNM clouds in the WNM, is £ P ~ 215 pc (§ 2). 

The parameter T = t C ooi/4hk that describes the strength of the turbulence is directly 
proportional to r tur b- Recall from § 2 that t coo i = |(l.ln/cT)/n 2 A for a neutral gas with 
10% He; in terms of the isothermal sound speed a t h, this is ^(fina^/ nA) . We estimated 
the shock time to be t s hk = ip/iy/S^th)- For a Kolmogorov-type spectrum (a 3 oc £), this 
becomes t sh k = (°"th/v / 2)(^/o' 3 )- We then find 



Thus, shocks are important in driving the gas away from thermal balance if and only if 
turbulent heating (evaluated with <p k = 1) is important. 

What do numerical simulations say about the value of the parameter $p. Stone et 
al. (1998) considered MHD turbulence driven by a range of wavelengths, with the power 
peaking at a length scale id- 5 Their results imply <f>^ = 0.94 for the case in which the initial 



5 In comparing with their work, it must be kept in mind that they expressed their results in terms of the 
size of the simulation box L, not £d — L/8, and they left open the question as to whether the dissipation 
rate is ~ 5v 3 /L or ~ Sv 3 /£d- We assume the latter, which is consistent with the subsequent work of Mac 
Low (1999). 
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field strength was such that the Alfven velocity va and the isothermal sound speed a t h were 
equal, as is approximately true in the diffuse ISM. In his study of turbulent dissipation, Mac 
Low (1999) calculated three models for the case va = Cth- He used two slightly different 
approaches for measuring the dissipation rate, one in terms of the volume-averaged rms 
velocity and one in terms of the mass-averaged rms velocity. Like Stone et al. (1998), we 
have used the latter approach. In our notation, Mac Low found 4>e = 0.9 ± 0.13 dex for 
this case, in fortuitously good agreement with Stone et al. (this agreement between the two 
calculations was not as good for other values of ct^/va)- 

In both these simulations, the turbulence is "balanced," in that the average wave power 
is the same in both directions along the field. Since the sources of interstellar turbulence are 
intermittent in space and time, actual interstellar turbulence is likely to be imbalanced, and 
Cho, Lazarian, & Vishniac (2002) have shown that this can substantially reduce the decay 
rate of the turbulence. The value 4>e ~ 1 found in these simulations is thus an upper limit 
to the value expected in the ISM. In our work, we shall somewhat arbitrarily adopt 0^ = 0.5 
as representative of IS turbulence. 

Recall that we parameterized the strength of the turbulence in terms of the ID veloc- 
ity dispersion at a scale of 1 pc, er(l). We estimated this from observation by assuming 
Kolmogorov scaling, which is reasonable since the motions of the CNM clouds observed by 
Heiles & Troland (2002) have a velocity dispersion that is very nearly the same as the ther- 
mal velocity in the WNM (for T = 8000 K, both are ~ 7.1 km s _1 ). In order to extend the 
estimate of the heating rate to other parts of the Galaxy, we assume that a is approximately 
constant and that t scales as the thickness of the H I disk, so that <r 3 /£ = cr(l; i?o) 3 oc 1/Hf 1 , 
where cr(l; Rq) is the value of cr(l) at the solar circle. From equation (B2) we find that the 
turbulent heating rate is then 

(B4) 

erg s" 1 H _1 (B5) 

Note that the estimated turbulent heating rate inside i? is constant, since the scale height 
of the H I does not change there. 
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B.2. Turbulent Energy from Differential Rotation 

Next, consider the extraction of turbulent energy from differential rotation. As pointed 
out by Sellwood & Balbus (1999), this process occurs at a rate 

( B r B <t>^ r \ dD. 
f*=-{— toT d^R' ( B6 ) 

where the rotation velocity is VLRcj) and the shear is in the radial direction. In a steady state, 
this rate of extraction of energy from differential rotation will be balanced by dissipation of 
energy, and so long as the velocities induced by the differential rotation are subsonic, much of 
this energy should be dissipated in a turbulent cascade. Hawley, Gammie, & Balbus (1995) 
have carried out MHD simulations of the generation of turbulent velocities and magnetic 
fields in a shearing box. They present detailed results for one model of a shearing box 
simulation; using these results, we find that the excitation rate simplifies to 

e = 0.755v 2 (—^-] • (B7) 

Note that this expression is based on the assumption that the velocities are generated by 
the differential rotation, as may be the case in the outer Galaxy; it may not apply if the 
velocities are generated by other mechanisms, such as supernovae. Even in regions where 
the turbulence is primarily due to differential rotation, this expression must be regarded as 
higly approximate since it does not take into account either the vertical structure or the 
multi-phase nature of the ISM. 

The dynamical model of the Galaxy that we have adopted from Dehnen & Binney 
(1998) has a rotational velocity that declines slowly beyond the solar circle, dQ/dlnR = 
(—244 km s~ 1 )/R for R > 8.5 kpc. We adopt 6 km s _1 as a typical velocity dispersion for 
gas in the outer parts of disk galaxies (e.g., Martin & Kennicutt 2001), somewhat less than 
the 7 km s _1 for the CNM and the 11 km s _1 for the WNM in the solar neighborhood (Heiles 
& Troland 2002). The turbulent heating rate due to differential rotation is then 

r turb = 1.50 x 10- 26 ( — °" ) 2 i- ergs" 1 !!- 1 (R k > 8.5 kpc). (B8) 
V 6 km s 1 / Rk 

This heating rate is within a factor 2 of that in equation (B5) for 8.5 < R k < 25 kpc. Given 
the uncertainties, these two estimates are in satisfactory agreement. The generic estimate in 
equation (B5) is about twice that for differential rotation in equation (B8) at the solar circle, 
consistent with the idea that supernovae are an important source of turbulent motions in 
the local ISM. 
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B.3. Results 

To estimate the effects of turbulent heating in the Galaxy, we adopt equation (B5), since 
it applies to both the case in which energy is injected by differential rotation and that in 
which it is injected by explosive events such as supernovae. It should be borne in mind that 
this estimate for the turbulent heating is quite uncertain, since it is based on highly idealized 
numerical simulations and an uncertain correction for an imbalanced turbulent cascade; on 
the other hand, it is reassuring that it is in accord with simple dimensional analysis (eq. B2 
with 0e of order unity). The estimate also depends on the amplitude of the turbulence in 
the subsonic regime, which is uncertain at present. Our estimate for the turbulent heating 
in the diffuse H I is about 30 times less than that invoked by Minter & Balser (1997) and 
Mathis (2000) in their studies of heating of the diffuse H II. 

Turbulent heating, unlike heating by cosmic rays, EUV/soft X-rays, and photoelectric 
heating, is independent of the depth into the cloud, and it therefore becomes more important 
at high column densities. If turbulent heating at the rate given by equation (B5) occurs 
throughout the Galactic disk, then for our fiducial column density of 10 19 cm -2 , at P m i n , our 
estimate of the turbulent heating rate exceeds the cosmic ray rate for Galactic radii between 
3 kpc and 18 kpc and amounts to 70% of the photoelectric heating rate at R k = 17 kpc. For 
iV = 10 20 cm -2 , the turbulent rate is always greater than the cosmic ray rate and is equal to 
93% of the photoelectric rate at R k = 17 pc. The average thermal pressure in a two-phase 
ISM including turbulent heating is given by 

Pth,avcA = 1.2 x 10 4 exp(-i4/7.5) K cm" 3 , (B9) 

where the fit is good to ±10% between 3 < Rk < 18 except at our R k = 11 model point 
where the fit overestimates the thermal pressure by 25% . 

The turbulent heating has the greatest effect in the outer Galaxy where heating rates 
based on stellar photons or supernovae are small. With turbulent heating, at R k = 17 kpc, 
we find P min /k = 650 K cm" 3 , P th ,aveA = 1370 K cm_3 > and Pmax A = 2900 K cm~ 3 which 
are factors of 1.7, 1.9, and 2.1 respectively times the non-turbulent case. We find that Pwnm' 
(the thermal pressure in the midplane based on the simple — and incorrect — assumption that 
the WNM is supported only by thermal pressure) falls below P max at R k = 15 kpc. The 
"mean" thermal pressure (Pwnm) (an estimate of the thermal pressure in the midplane based 
on the assumptions that all the H I is WNM and that the WNM fills space) falls below P max 
at approximately R k = 13.5. It follows that the H I cannot be all WNM out to R k = 13.5, 
and a two-phase medium must exist out to that point. Since (Pwnm) remains above P min 
out to R = 18 kpc, we conclude that pressure fluctuations will in fact produce a two-phase 
medium to these distances. 
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C. Analytic Thermal Balance Model for Cool Gas 

When gas is in thermal balance, the thermal pressure P can be expressed as a function 
of either density or temperature. We seek P(T) for the temperature range T ^ 1000 K where 
[C II] 158 /iin and [O I] 63 /im radiation dominates the gas cooling, and for n < n^ r ([C II]) ~ 
3000 cur 3 and n e < n^ r ([C II]) ~ 30 cur 3 so that the [C II] + [O I] cooling is proportional 
to n 2 . 



C.l. Thermal Balance: Heating and Cooling 

Cooling rates per unit volume can be written as n c nA c s1 where n c is the density of the 
collisional agent (n e or n is the hydrogen nucleus density, and A^ is the cooling rate 
coefficient, which takes into account the gas phase abundance of the species at the solar 
circle (Table 1). We find 

Ag „ = 3.15 x 10- 27 e~°- 92 ^Z' g erg cm 3 s" 1 , (CI) 

A e c n = 1.4 x 10- 24 T 2 - 1/2 e- a92 / T2 4 erg cm 3 s" 1 , (C2) 

and 

Ag ! = 2.5 x 10- 27 T^ A e- 2 - 28 ^Z' g erg cm 3 s" 1 , (C3) 

where T2 = T/(100 K). We assume that H atom collisions dominate electron collisions. 
Comparison of uhAq n with n e A c n shows that this condition is equivalent to x e = n e /n < 
2.3 x 10- 3 T 2 1/2 . 

Assuming that H atoms dominate the excitation of [C II] and [O I], that n < n cr ~ 3000 
cm~ 3 , and 100 K < T < 1000 K (the relevant regime for T( min )), we use numerical results 
using the above cooling coefficients Ag n and Aq i to obtain a simple form for the total 
cooling coefficient 

Af ot = 5.4 x 10- 27 T 2 a V L5 / T2 ^ erg cm 3 s^ 1 , (C4) 

which is good to ±20%. A simple form such as this is required, rather than Aq u + Aq l , in 
order to analytically determine T( min ). 

We modify the grain photoelectric heating rate of Bakes & Tielens (1994) by multiplying 
by a factor of 1.3, which accounts for a higher PAH abundance (6 x 10~ 7 by number relative 
to hydrogen) compared with that assumed in Bakes & Tielens (1994). The heating rate per 
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unit volume for T ^ 1000 K is given approximately [see § 5, eq. (23)] 

1.1 x lO-^G'nZ'n , , 

nr pc = ° * . 73 erg cnT 3 s" 1 . (C5) 

1 + 3 2 x 10~ 2 I G '° T ' 2 I 

l n c (f> PAH J 

Recall (see § 5) that <Ppah is a parameter of order unity which scales the PAH collision rates; 
4>pah = 0.5 in our standard model. 

Note that when the unity term dominates in the denominator, the photoelectric heating 
is not significantly suppressed by positive charging, i.e., the second term in the denominator 
represents effects of grain charging. To obtain an analytic solution for T( min ) below, we will 
require the second term to dominate, as it often does for a range of conditions centered on 
solar neighborhood values. 

C.2. The Electron Density n e 

In order to obtain an analytic expression for P under thermal balance conditions, an 
analytic expression for n e is required to substitute into the grain photoelectric heating equa- 
tion (C5). The chemistry leading to steady state electron abundances is quite complex and 
interesting. We use the results of the numerical code to determine the dominant reaction 
chains leading to n e , and to determine the major competitors to these reactions. We then 
find an analytic expression for n e , along with the conditions required to ensure the assumed 
reaction chain. 

For densities and temperatures near n min and T( min ) at the solar circle, H + is the domi- 
nant ion species and therefore n c ~ n H +- H + is produced mostly by the EUV/soft X-ray pho- 
toionization of H with a minor contribution from cosmic rays. We find that khc+ ~ 0.3nn+ 
under a variety of conditions at a column of N c \ ~ 10 19 cm -2 , a result due to the higher pho- 
toionization cross section of He at soft X-ray energies counterbalanced by the 0.1 abundance 
of He relative to H. Since H + is a surrogate to obtain the electron abundances, we roughly 
account for He + by increasing the photoionization rate of H by 30% to a rate 1.3 x 10~ 16 
s _1 at N c \ = 10 19 cm -2 appropriate to the solar neighborhood. The rate Q is expressed as 
Q = ( t /10~ 16 s _1 , where ( t is the total ionization rate (including primary and secondary ion- 
izations) of H by photons and cosmic rays. The destruction of H + is dominated by reactions 
with PAH~ or 

H + + PAH" -> H + PAH , «i = 8.3 x l()- 7 (f) PAH T 2 - 5 cm 3 s" 1 , (C6) 
where K\ is the rate coefficient calculated from Draine & Sutin (1987) using disk PAHs with 
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N c = 25 carbon atoms and the disk radius a = (iVc/1.222) - 5 . Competing reactions are 

H + + PAH -> PAH+ + H , k 2 = 3.1x lO~ 8 PA //T 2 - a5 cm 3 s" 1 , (C7) 

and 

H+ + e -> H + hv , k 3 = 8.0x lO" 12 ^" - 75 cm 3 s" 1 , (C8) 

where the reaction coefficients k 2 and k 3 are calculated similarly to K\. Equating the forma- 
tion of H + to the destruction of H + , we obtain for the electron density 

1.3xlO- 16 C' n , „ , CTl /2 n 3 
n e = ^— = 1.6 x l(r 10 --^ cm" 3 . (C9) 

^l^PAH- 0PAH^PAH- 

The analytic solution for n c therefore requires a solution for u P ah- ■ We simplify PAH 
ionization state calculations by assuming single-sized PAHs with three ionization states, 
PAH~, PAH , and PAH + . We find numerically that, over a wide range of relevant parameter 
space, neutral PAH dominates the population with tipah — 0.7npAH where npAH — 6.0 x 
10~ 7 nZ' d is the total density of PAHs in all states. PAH + is formed by the FUV photoreaction 

hv + PAH -> PAH+ + e , k 4 = 7.85 x lO" 9 ^ s" 1 , (CIO) 

where the rate coefficient is calculated from Bakes & Tielens (1994) using Nq = 35 carbon 
atoms and assuming a disk geometry. PAH + is destroyed primarily by recombination with 
electrons 

PAH+ + e -> PAH , k 5 = 3.5 x lO^p^T^ - 5 cm 3 s" 1 . (Cll) 
The condition that npAH+/^PAH° < 1 is then 7.85 x 10~ 9 G' /(K, 5 n e ) < 1 . PAH~ is formed by 

PAH + e -> PAH~ , k 6 = 1.34 x 10- 6 <f> PAH cm 3 s" 1 , (C12) 

and PAH~ is destroyed primarily by FUV photodetachment 

PAH" + hv -> PAH + e , k 7 = 2.00 x lO -8 ^ s" 1 . (C13) 

The condition n P AH-/^PAH () < 1 is then /t 6 n e /2.0 x 10~ 8 Gq < I. 

Assuming npAH = 0.7npAH, we find (with n and n c in units of cm~ 3 ) 

npAH- = 2.8 x \§- h nMpAHZ' A G'^ cm" 3 , (C14) 

or substituting into equation C9, we find 

n e = 2.4 x 10- 3 C t /1/2 T 2 1/4 ^/ 2 Z7- 1/2 ^ cm" 3 . (CIS) 
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C.3. Analytic Expression for P(T) and Conditions for Validity 

We substitute equation (C15) for n e into equation (C5) for the heating, equate the 
heating to the cooling given by equation (C4), and obtain for neutral gas in thermal balance 
with 100 K < T < 1000 K 

2240 G' ( §) T 2 - 8 e L5 /^ 

P/k = l.lnT = ) s \ /2 x 0.365 • (C16) 

1 + 2.6 ° % d 



This equation is valid if all of the following conditions are met: 

(1) ^° 2 d < 17 (C17) 

Ct 

Condition (1) ensures that H + is mostly destroyed by PAH~, and not PAH . This condition 
and the rest below are derived by using the rate coefficients above, along with equation 
(C15) for n e . 

(2) U ° l2 Zd < 112 (C18) 

Q 

Condition (2) ensures that neutral PAH and not PAH + dominates the PAH population. 
Note that if condition (1) is satisfied, condition (2) is automatically satisfied. 



(3) ' 0.!>,A>7.. W/ (C19) 



Condition (3) ensures that H atoms, and not electrons, dominate collisional excitation of 
[C II] 158 fim. In these equations, the hydrogen nucleus number density n is in cm -3 . 

(4) < 2.9ru^ PAH (C20) 
Condition (4) ensures that H + is destroyed by PAH - , and not by electron recombination. 

(5) ^± > 2.5 x 10- 2 (C21) 

t 2 1/2 c; 

Condition (5) ensures that neutral PAH , and not PAH" dominate the PAH population. 

(6) G '° Z 'f 2 > 9.5 x lO" 3 (C22) 
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Condition (6) ensures that FUV dominates the destruction of PAH , and not reactions with 
H + . Note that for T 2 ~ 1 — 10, condition (5) automatically insures condition (6). 

fll f-lrp^-l'i 

(7) > 3.5 x (C23) 

Condition (7) ensures that the dominant ion is H + and not C + , which is ignored in the 
analytic analysis. These conditions ensure the validity of equation (C16) for P(T). 

Although the number of conditions is large, inspection of each of them reveals that for 
n ~ 10 cm -3 and <Ppah — 0.5, all of the conditions are met for a wide range of conditions 
centered on solar neighborhood values. The motivation for deriving P(T) is to find an 
analytic solution for P m i n , TJmin), and n min . Thus, the range of validity is centered on T ~ 
T( min ) ~ 250 K and n ~ n min ~ 10 cm -3 , as we shall derive below. To calculate P m i n , we 
take dP/dT = and solve for T( min ). Substitution of T( min ) into equation (C16) gives P m i n . 
No analytic solution is possible unless we simplify the denominator of P(T) (eq. [C16]). 
If conditions (1-7) are satisfied, then the second term in the denominator is usually larger 
than unity. The denominator originates from the grain photoelectric heating equation, and 
the second term corresponds to the effects of positive charging on the grains/PAHs. The 
condition that the second term dominates is 

ft l yi rpl/2 

(8) d , 2 > 0.072 . (C24) 

Ct 

Note that condition (8) and condition (5) are nearly identical, for T 2 ~ 2.5. With all eight 
conditions satisfied, we find the solutions for T( min ), P m i n , and n min given in the main text. 
Substitution of n min (with the factor of unity in the denominator removed) and T( min ) into 
the eight conditions give the simplified set of conditions given in the text. 
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Table 1. Model Parameters - Local Values a 



Parameter Value 



Reference 



J FUV ( J R ) b 2.2 x 10" 4 erg cm~ 

Ccr c 1.8 x 1(T 17 s- 1 

CxR d 1.6 x 10- 17 s- 1 

A c e 1.4 x 10~ 4 

A) f 3.2 x 10~ 4 

^ PAH g 6 x 10~ 7 

A y h N/(2 x 10 21 cm" 2 ) 

Shi 1 5 M pc~ 2 

S H2 j 1.4 M pc- 2 

S W NM k 2.75 M pc~ 2 

(nni) 1 0.57 cm" 3 

i^l 1 Im 115 pc 

H^ n 59 pc 



i 



sr 



-i 



5,6 
7 

4,8 
9 

10 
11 
10 
10 
10 

11 



1,2 



3 

4 



Refs.-(l) Draine (1978), (2) Habing (1968); (3) WHMTB; (4) This paper; (5) Cardelli et al. (1996); 
(6) Sofia et al. (1997); (7) Meyer et al. (1998); (8) Tielens et al. (1999); (9) Bohlin, Savage, & Drake 
(1978); (10) Dickey & Lockman (1990); (11) Bronfman et al. (2000) 
a Values at solar circle R = Rq = 8.5 kpc. 

b Intensity of FUV (6 eV < hv < 13.6 eV) interstellar radiation field. This intensity is a factor 1.7 
higher than the integrated field of Habing (1968) which has a value of J FUV (i?o) = 1.3 x 10 -4 erg 
cm~ 2 s _1 sr^ 1 . The interstellar field has a value of Go = 1.7 in units of the Habing field. 
c Primary cosmic-ray ionization rate of hydrogen. 

d Primary EUV plus soft X-ray ionization rate of hydrogen at a cloud depth of iV c i = 10 19 cm~ 2 . 
The primary rate depends on the adopted value of the WNM cloud column, N c \. The total rate 
is higher than the primary rate due to secondary ionizations and increases with lower electron 
fraction. For our standard model, over the 2 phase range, the total rate is a factor 3 — 5 higher 
than the primary rate. For iV cl = 1.0 x 10 20 cm" 2 , Cxr(^o) = 8.9 x 10~ 19 s _1 ; N d = 3.0 x 10 18 
cm" 2 , Cxr(^o) = 7.6 x lO" 17 s~\ 
e Gas phase carbon abundance per H nucleus. 
f Gas phase oxygen abundance per H nucleus. 

g PAH abundance per H nucleus. This abundance gives a total number of C atoms in PAHs of 
22 x 10 -6 relative to hydrogen. 

h Magnitudes of visual extinction per hydrogen column density. 

Atomic hydrogen surface density through full disk, iV(H I) = 6.25 x 10 20 cm -2 . 

^Molecular hydrogen surface density through full disk, 2iV(H2) = 1.75 x 10 20 cm~ 2 . 

k WNM column density through full disk, iVwNM,d = 3.45 x 10 20 cm~ 2 . 

l Mean H I density. 



-73- 



111 H I half height to half intensity. 
n H2 half height to half intensity. 



Table 2. Model Parameters - Galactic Values 3, 



R 


G' h 


(-1 c 
SCR 


{i d 

SXR 


Z' A = z' g ° 


(kpc) 










3 


2.95 


3.64 


5.02 


2.43 


4 


2.73 


2.56 


3.62 


2.07 


5 


2.21 


1.83 


2.72 


1.76 


8.5 (= Ro) 


1.00 


1.00 


1.00 


1.00 


11 


0.509 


0.424 


0.359 


0.668 


15 


0.198 


0.194 


0.156 


0.351 


17 


0.116 


0.182 


0.145 


0.254 


18 


0.0634 


0.176 


0.140 


0.216 



a Values scaled to solar circle. 

b Scaled Intensity of FUV interstellar radiation field. 

c Scaled primary cosmic ray ionization rate 

d Scaled primary EUV and soft X-ray ionization rate. 

e Scaled gas phase metallicity, Z' , and dust/PAH abundances, Z' d 



Table 3. Range of Physical Conditions for Two-Phase Medium 
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Table 4. Dependence on Model Parameters* 



Model 


Standard 13 


Low 0pah c 


High 0pah d 


Low np\n/n e 


Low Gq 


P min (K cm- 3 ), N = N c f 


1960 


1560 


2270 


1580 


1460 


P max (K cm" 3 ), N = iV c i g 


4810 


3150 


5970 


3920 


3980 


P th , avc h (K cm- 3 ), iV = 7V cl s 


3070 


2220 


3680 


2490 


2410 


T ave i (K),iV = AT cl s 


85 


94 


96 


79 


76 


n avc J(K), N = N cl z 


33 


21 


35 


29 


29 


nA k (erg s- 1 H- 1 ), N = JV d « 


4.1(-26) 


3.2(-26) 


5.0(-26) 


3.3(-26) 


3.1(-26) 


C F/C Itot 1 , iV = iV cl s 


0.190 


0.145 


0.205 


0.168 


0.167 


T^ 1 (K), 1 x 10 20 (cm" 2 ) 1 " 


71 


67 


84 


64 


65 


nA k (erg s- 1 H- 1 ), 1 x 10 20 (cm" 2 ; 


) m 3.6(-26) 


2.6(-26) 


4.3(-26) 


2.8(-26) 


2.7(-26) 


C I*/C Itot 1 , 1 x 10 20 (cm" 2 ) 111 


0.201 


0.162 


0.213 


0.178 


0.172 


C I/C IP , 1 x 10 20 (cm- 2 ) m 


2. 3 (-3) 


7.1(-4) 


6.1(-3) 


1.5(-3) 


3.7(-3) 



a a(— b) means a x 10" 6 . 

b Standard model with </> PA h = 0.5, n PAH /n = 6 x 10" 7 , G = 1.7, N d = 1 x 10 19 cm" 2 . 
c Low PAH collision rates with <^> PAH = 0.25, n PAH /n = 6 x 10" 7 , G = 1.7, JV d = lx 10 19 cm" 2 . 
d High PAH collision rates with (p PA H = 1.0, np A n/n = 6 x 10" 7 , G = 1.7, N d = 1 x 10 19 cm" 2 . 
e Low PAH abundance with 0pah = 0.5, np A n/n = 4 x 10" 7 , Go = 1-7, A^i = 1 x 10 19 cm" 2 . 
f Low FUV field with PA h = 0.5, n VAn /n = 6 x 10" 7 , G = 1.1, N cl = 1 x 10 19 cm" 2 . 
g Model result at cloud depth of A^i = 1 x 10 19 cm" 2 and Pth,ave- 

Pth, ave — \/Pmax X P mm . 

'Temperature of CNM at Pth, ave- 
j Density of CNM at P th ,ave- 

k Gas cooling rate per hydrogen atom from [C II] 158 fj,m line emission in CNM. 
X C I*/C Lot population ratio. 

111 Model result at cloud interior at depth of 1 x 10 20 cm" 2 and Pth, ave- 
n C I/C II abundance ratio. 
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Table 5. Parameters in 


Galactic 


Disk Potential Moc 


lel a 


Component 


2, 


Rm 


Rd 


z d 




(M pc- 


~ 2 ) h (kpc) 


(kpc) 


(pc) 


H I (R < 13 kpc) 


7.94 


1.0 


1000 


178 


H I (R > 13 kpc) 


571 


10 


4.00 


324 


H 2 


57.5 


3.3 


2.89 


63.4 


H II 


1.39 





30.0 


880 


Thin star disk {R < R ) 


1058 





2.55 


180 


Thin star disk (R > R ) 


1058 





2.55 


z d 


Thick star disk 


70.6 





2.55 


1000 



a Disk potential model based on Dehnen & Binney (1998). Each component has a density- 
distribution given by the form p(R) = T ld (2z d )^ 1 exp(— [R rn /R] — [R/R d ] — [\z\/z d ]). Bulge 
and halo components are the same as Dehnen & Binney (1998) model 2b. 
b Includes He mass. 

c Stellar disk height given by z* d = 180 x E(i? )/E(i?) pc. 



-77- 



Figures 

Fig. 1. — Azimuthally averaged H I surface density Shi (dash) and H 2 surface density £h 2 
(dot) in the Galactic disk versus Galactocentric radius R. Total H I plus H 2 is shown as a 
solid line. Mass does not include He. 

Fig. 2. — Azimuthally averaged H I half width to half maximum height H^ 1 versus Galac- 
tocentric radius R. 

Fig. 3. — Mean H I density in the Galactic midplane versus Galactocentric radius R. The 
value at R = R is taken to be (nui) = 0.57 cm" 3 and scaled by (n H i) oc Y.m/Hf 1 at other 
radii. 

Fig. 4. — Calculated FUV opacity in the Galactic midplane versus Galactocentric radius R. 

Fig. 5. — Calculated FUV field in the Galactic midplane versus Galactocentric radius R 
normalized to the value at R = R . At R = Rq the field strength is equal to 47rJ FUV (i? ) = 
2.7 x 10~ 3 ergs cm" 2 s _1 . 

Fig. 6. — Primary cosmic-ray (solid) and EUV plus X-ray (dot) ionization rates versus 
Galactocentric radius R normalized to the value at R = Rq. At R = Rq the primary cosmic- 
ray ionization rate is taken to be Ccr(-Ro) — 1-8 x 10~ 17 cm" 3 s" 1 . The EUV and X-ray rate 
depends on the adopted value of the WNM cloud column, N c \. For N c \ = 1.0 x 10 20 cm" 2 , 
Cxr(-Ro) = 8.9 x 10" 19 s" 1 ; N cl = 1.0 x 10 19 cm" 2 , ( X r(R ) = 1-6 x 10" 17 s" 1 ; N cl = 3.0 x 10 18 
cm" 2 , Cxr(-Ro) = 7.6 x 10" 17 s" 1 . Note that the EUV and X-ray rate always exceeds the 
cosmic ray rate because of the effects of secondary ionizations. Typical values for N c \ are of 
order 10 19 cm -2 . 

Fig. 7. — Phase diagrams showing thermal pressure Pj k versus hydrogen nucleus density n 
at Galactocentric radii R = 3,5, 8.5, 11, 15, and 18 kpc. Curves apply to the WNM/CNM 
boundary at a depth of 1 x 10 19 cm" 2 through the WNM. Gas is thermally stable to isobaric 
perturbations where dP/dn > 0. 

Fig. 8. — Thermal pressure P/k versus hydrogen nucleus density n at R — 8.5 kpc showing 
the effects of varying the collision rate parameter 0pah and the PAH abundances. For the 
standard PAH abundances (np^/n = 6 x 10" 7 ), curves are shown for 0pah — 0.25 (dash), 
0pah = 0.5 (solid), and 0pah = 1-0 (dot). Our standard model uses 0pah = 0.5. Also shown 
is a curve for 0pah = 0.5 and a low PAH abundance of npAii/n = 4 x 10" 7 (long dash). 

Fig. 9. — Thermal pressure P/k versus hydrogen nucleus density n at Galactocentric radius 
R = 8.5 kpc. Curves are shown for various values of WNM atomic column density N c \ and 
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apply to the WNM/CNM boundary. 

Fig. 10. — Upper Panels: Heating and Cooling curves versus hydrogen nucleus density n 
at various Galactic distances, R. Heating rates (dash); Photoelectric heating from small 
grains and PAHs (PE); EUV and X ray (XR); Cosmic ray (CR); photoionization of C (C I). 
Cooling rates (solid); C II 158 jum fine-structure (C II); O I 63 jum fine-structure (O I); 
Recombination onto small grains and PAHs (Rec); Lya plus metastable transitions (Lya); 
C I fine-structure 609 fxm (C I*); C I fine-structure 370 /im (C I**). Lower Panels: Gas 
temperature T (solid) and electron fraction n e /n (dash) versus hydrogen nucleus density n. 
(a) R = 5 kpc. (b) R = 8.5 kpc. (c) R = ll kpc. (d) R = 17 kpc. 

Fig. 11. — Infrared intensity in the Galactic midplane versus Galactic longitude. COBE 
observations (solid) from Sodroski et al. (1994). Calculated dust emission (dot) using the 
FUV opacity from Fig. 4, and the FUV field from Fig. 5. 

Fig. 12. — Limiting thermal pressures in the Galactic midplane versus Galactocentric radius 
R. -Pwnm' (solid) is the thermal pressure produced by the H I layer assuming all of the 
gas is in the form of WNM and supported in hydrostatic equilibrium by thermal pressure. 
-Pmax (short dash) and P m i n (long dash) are the maximum and minimum pressure range for 
a two-phase medium. Pth,avc (dash) is the mean pressure Pth.ave = (-Pmin x Pmax) 1 ^ 2 - Panels 
show the effects of different WNM columns, (a) N d = 1.0 x 10 19 cm' 2 , (b) N cl = 1.0 x 10 20 
cm" 2 , (c) N d = 3.0 x 10 18 cm" 2 , (d) N d = 1.0 x 10 19 cm" 2 curves with (P W nm) added. 
(Pwnm) is the thermal pressure that would be present in the Galactic midplane if all of the 
H I had a temperature of 8000 K and a density given by the mean density shown in Fig. 3 
(see § 6.4.2). 

Fig. 13. — Region of validity for the analytic solution of P m i n (eq. [37]) as a function of the 
dust abundance Z' d , gas phase metal abundance Z' g , ionization rate Q, and FUV radiation 
field G' . All parameters are scaled to their value in the solar neighborhood. The valid region 
is shown for the case of dust and metal abundances scaling linearly with the elemental 
abundance (Z' = Z' d = Z'). The shaded region shows the range in which the analytic 
solution is good to within ±50% of the numerical results. Also shown are the more restrictive 
conditions derived in Appendix C and given in equations (40) through (42): (solid) PAHs are 
positively charged and the destruction of H + is dominated by reactions with PAH" rather 
than with PAH (eq. [40]); (dot) Electrons are supplied by H + rather than C + (eq. [41]); 
(dash) the destruction of H + is dominated by reactions with PAH" rather than by electron 
recombination (eq. [42]). The analytic solution is good to within ±45% of the numerical 
results within this region. 
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